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SUMMARY 


The  technique  of  Large-Eddy  Simulation  (LES)  is  applied  to  flow  over  terrain,  in 
order  to  generate  turbulent  velocity  fields  for  use  in  dispersion  calculations.  The 
objective  of  the  research  is  to  study  the  effects  of  simple  terrain  on  the  diffusion  of  a 
tracer  from  a  localized  source.  Most  of  our  understanding  of  atmospheric  turbulence  is 
based  on  studies  over  horizontally-homogeneous  terrain,  but  the  presence  of  hills  or 
ridges  disturbs  the  atmospheric  boundary  layer  and  presumably  affects  the  dispersion 
properties.  LES  provides  the  most  fundamental  methodology  available  for  studying 
atmospheric  turbulence,  and  is  therefore  particularly  appropriate  for  modeling  dispersion 
in  the  presence  of  idealized  terrain. 

We  first  evaluate  the  capability  of  LES  using  laboratory  data  on  neutral  flow  over 
sinusoidal  ridge  terrain  with  a  variety  of  amplitudes.  Good  agreement  is  obtained  for 
both  mean  flow  and  observed  turbulence  statistics.  We  then  utilize  the  LES  model  to 
study  atmospheric  boundary  layer  flow  over  similar  sinusoidal  terrain,  and  calculate  the 
dispersion  of  a  tracer  using  a  Lagrangian  puff  model  and  also  a  particle  trajectory  model. 
The  effects  of  heat  flux  on  the  dispersion  properties  are  examined,  and  found  to  reduce 
the  effect  of  the  terrain.  This  is  due  to  the  increased  vertical  mixing  under  convective 
conditions,  which  limits  the  effect  of  the  increased  shear  near  the  terrain  surface. 

Quantitative  information  on  dispersion  rates  has  been  obtained  for  idealized 
terrain  profiles.  Several  basic  dispersion  phenomena  associated  with  flow  and  turbulence 
fields  over  terrain  have  been  elucidated.  In  addition  to  the  effects  of  convection  on 
streamwise  dispersion,  the  LES  calculations  also  show  the  effect  of  separated  flow 
regions.  Material  released  or  entrained  into  a  separated  flow  region,  such  as  that  within  a 
steep-sided  valley,  is  detrained  over  a  relatively  long  timescale.  The  detrainment  is  not 
incremental,  however,  but  intermittent  in  nature  on  the  timescale  of  the  recirculation. 
This  implies  that  a  small  release  within  a  separated  flow  will  be  randomly  ejected,  and 
the  downwind  concentration  is  therefore  extremely  uncertain. 

In  summary,  LES  has  been  shown  to  be  a  valuable  tool  for  studying  boundary 
layer  flow  and  dispersion  over  simple  terrain. 
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1.  INTRODUCTION 


The  presence  of  terrain  can  greatly  influence  the  structure  of  the  atmospheric  boundary 
layer,  modifying  both  the  mean  flow  and  the  turbulent  fluctuations  from  the  values  over  a  flat 
surface.  Flow  acceleration  over  hillcrests,  separation  and  flow  reversal  behind  steep  hills,  and 
buoyancy-driven  drainage  circulations  are  all  familiar  examples  of  flow  phenomena  associated 
with  topography.  Many  studies  have  been  conducted  to  investigate  turbulent  boundary  layer 
flow  over  terrain,  and  the  wide  range  of  phenomena  has  necessitated  many  simplifying 
assumptions.  Since  the  complexity  of  the  real  atmospheric  situation  masks  the  individual 
processes,  most  theoretical  studies  focus  on  idealized  terrain  to  isolate  specific  aspects  of  the 
dynamics.  Thus,  for  example,  many  researchers  have  studied  neutral  or  near-neutral  flow  over 
topography,  using  a  variety  of  turbulence  closure  models,  to  elucidate  the  effect  of  different 
empirical  modeling  assumptions  on  the  predictions.  For  example,  Jackson  and  Hunt  (1975), 
Mason  and  Sykes  (1979),  and  Hunt  et  al.  (1988)  present  linearized  analytical  solutions  using  a 
turbulent  eddy-viscosity  model  for  flow  over  idealized  terrain.  Higher-order  turbulence  closure 
models  have  been  applied  by  Sykes  (1980),  Newley  (1986),  and  Belcher  et  al.  (1993),  and 
numerical  solutions  for  various  closure  models  have  been  presented  by  Taylor  et  al.  (1976)  and 
Wood  and  Mason  (1993).  Belcher  and  Hunt  (1998)  review  the  subject  with  emphasis  on  the 
turbulence  distortions  and  the  interaction  with  vertical  wind  shear. 

We  will  focus  on  the  aerodynamic  flow  studies  described  above,  since  they  are  concerned 
with  smaller  scale  terrain  features  and  their  effects  on  boundary  layer  turbulence.  This  is  directly 
relevant  to  the  dispersion  of  material  from  localized  sources,  and  extends  the  homogeneous 
terrain  approximations  used  in  traditional  models  for  the  turbulent  boundary  layer  structure. 
Other  researchers  have  examined  mean  flows  induced  by  larger  scale  terrain  features  (mountain 
ranges)  and  buoyancy-driven  circulations  (valley/slope  winds),  but  these  studies  are  principally 
concerned  with  mean  flow  effects,  and  place  less  emphasis  on  the  turbulence  modifications. 
Since  our  main  interest  is  in  extending  the  atmospheric  diffusion  parameterizations  to  include 
small-scale  terrain  effects,  we  will  concentrate  on  terrain  with  length  scales  comparable  to  the 
atmospheric  boundary  layer  depth,  i.e.,  on  the  order  of  a  few  kilometers  or  less. 

The  terrain-induced  distortions  of  the  flow  imply  modification  of  the  dispersion 
properties  of  the  boundary  layer,  since  material  is  transported  by  the  turbulent  wind  field.  Here, 
it  is  important  to  distinguish  mean  flow  effects  from  turbulent  diffusion  effects,  although  the  two 
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will  be  coupled  in  general.  The  dispersion  properties  of  the  modified  boundary  layer  turbulence 
have  not  been  extensively  studied,  and  will  therefore  be  the  major  objective  of  the  current  work. 

As  discussed  above,  there  have  been  a  number  of  approaches  using  turbulence  closure 
models  to  calculate  the  effect  of  terrain  on  the  turbulent  boundary  layer.  Large-eddy  simulation 
(LES)  is  an  alternative  approach  to  Reynolds-average  closure,  which  minimizes  the  empirical 
assumptions  involved  in  the  calculation.  LES  uses  a  3-dimensional  time-dependent  numerical 
solution  of  the  spatially  averaged  Navier-Stokes  equations  to  represent  the  larger  scale  turbulent 
eddies  explicitly.  Turbulent  motions  smaller  than  the  resolved  scale  must  still  be  parameterized 
using  some  form  of  closure  model,  since  it  is  impractical  to  resolve  the  turbulent  spectrum  down 
to  viscous  scales  in  the  atmosphere.  However,  predictions  are  less  sensitive  to  the  subgrid 
closure,  and  the  dependence  is  reduced  as  the  numerical  resolution  is  improved. 

Recent  studies  of  turbulent  boundary  layer  flow  over  terrain  have  been  reported  by 
Krettenauer  and  Schumann  (1992),  for  convective  conditions,  and  Gong  et  al.  (1996),  for  neutral 
conditions.  Both  papers  study  two-dimensional  sinusoidal  terrain,  since  this  is  the  simplest 
possible  geometry  that  can  be  considered  using  the  periodic  lateral  boundary  conditions  used  in 
most  LES  calculations.  Krettenauer  and  Schumann  consider  zero  mean  wind  and  focus  on  the 
effects  of  terrain  on  the  convective  eddies,  while  Gong  et  al.  study  neutral  flow  over  waves  and 
compare  with  wind  tunnel  observations.  In  this  report,  we  will  examine  wind  flow  over  simple 
terrain  with  specific  emphasis  on  the  dispersion  of  a  passive  tracer  in  the  flow.  However,  we 
first  compare  the  LES  results  with  available  laboratory  data  on  flow  over  sinusoidal  terrain,  in 
order  to  examine  the  ability  of  LES  to  model  the  terrain-induced  phenomena. 
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2.  COMPARISON  WITH  EXPERIMENTAL  DATA 


Since  Large-eddy  simulation  of  flow  over  terrain  is  not  well-established,  we  first  make  an 
extensive  comparison  between  detailed  laboratory  measurements  and  the  numerical  predictions. 
Prior  to  examining  dispersion  of  a  passive  tracer  within  the  numerically  generated  turbulent  flow 
field,  we  must  ascertain  that  the  computational  results  are  consistent  with  real  flows.  To  this 
end,  we  have  compared  LES  predictions  with  a  number  of  laboratory  experiments,  as  reported  in 
Henn  and  Sykes  (1999). 

An  extensive  set  of  measurements  of  channel  flow  over  wavy  surfaces  has  been  presented 
by  Hanratty  and  coworkers  at  the  University  of  Illinois.  These  experiments  cover  a  wide  range 
of  slopes  and  Reynolds  numbers  since  the  nature  of  the  flow  depends  strongly  on  these 
parameters.  We  have  concentrated  on  the  experiments  with  larger  Reynolds  numbers  and  steeper 
slopes  since  they  are  most  relevant  to  the  atmospheric  flows  we  are  interested  in.  However,  we 
have  made  some  comparisons  over  a  range  of  slopes,  including  a  flat  channel. 

In  modeling  the  laboratory  and  atmospheric  flows,  the  major  required  modification  to  the 
basic  Cartesian  LES  technique  presented  previously  in  Sykes  and  Henn  (1989,  1992)  involves 
the  use  of  the  terrain-following  coordinate  transformation  of  Clark  (1977).  This  introduces  a 
number  of  additional  terms  to  the  basic  equations  and  complicates  the  numerical  implementation. 
In  particular,  the  pressure  field  must  be  computed  iteratively,  with  the  number  of  iterations 
increasing  with  slope.  Thus  the  computational  effort  is  greatly  increased.  Details  can  be  found 
in  Henn  and  Sykes  (1999). 

Although  the  laboratory  experiments  involve  fully  turbulent  flow,  finite-Reynolds- 
number  effects,  especially  near  the  lower  wavy  boundary,  are  still  evident.  Therefore,  the  basic 
LES  model  is  modified  for  this  comparison  with  the  implementation  of  a  no-slip  boundary 
condition  and  a  wall-damping  function  on  the  turbulent  length  scale.  Thus,  the  wall  shear-stress 
is  determined  by  the  viscosity  and  the  assumption  of  parallel  flow  at  the  first  grid  level.  The 
momentum  and  turbulent  kinetic  energy  transport  equations  are  also  modified  with  viscous 
diffusion  terms,  but  it  should  be  emphasized  that  it  is  not  intended  to  construct  a  true  low- 
Reynolds-number  model  (e.g.,  one  that  would  handle  transitional  flows),  but  is  included 
principally  to  impose  the  no-slip  boundary  condition. 
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This  study  involves  model  comparisons  with  experiments  over  a  wide  range  of  wave 
slopes,  from  flat  to  moderate  (with  nearly  linear  perturbations)  to  steep  (with  strongly  separated 
flow).  The  available  measurements  include  mean  streamwise  velocity  and  its  variance,  skewness 
and  flatness,  wall  shear-stress,  surface  pressure  and  integrated  form  drag.  Although  the 
calculations  use  moderate  resolution  (typically  48  points  in  each  horizontal  direction  and  around 
70  in  the  vertical),  the  LES  model  appears  capable  of  predicting  most  of  the  salient  features  of 
the  laboratory  flow.  In  particular,  the  response  of  the  mean  flow  and  streamwise  velocity 
variance  to  different  wave  amplitudes  is  accurately  predicted  even  close  to  the  surface.  Model 
sensitivity  to  resolution  and  domain  size  confirms  the  moderate  resolution/domain  results.  The 
resolution  tests  also  indicate  that  the  sub  grid  model  results  in  the  appropriate  transfer  of  energy 
to  the  resolved  scale  as  resolution  increases. 

One  striking  aspect  of  the  LES  results  is  a  large  increase  in  lateral  velocity  fluctuation 
intensity  on  the  upslope  of  the  ridges.  The  increase  seems  to  be  associated  with  vortex  stretching 
as  the  flow  is  accelerated  over  the  crest  of  the  ridge,  and  the  vortices  are  relatively  long-lived. 
The  increase  in  the  turbulence  intensity  scales  with  the  square  of  the  terrain  slope,  indicating  that 
it  is  a  nonlinear  phenomenon  in  contrast  to  the  linear  distortions  of  the  “upstream”  turbulence. 
Unfortunately,  the  lateral  velocity  component  was  not  measured  in  any  of  the  laboratory 
experiments,  so  this  feature  is  still  to  be  confirmed. 

For  applications  to  atmospheric  flows,  it  is  important  to  be  able  to  predict  variations  in 
the  mean  velocity,  turbulence  and  surface  pressure  along  the  wave,  as  well  as  to  characterize 
separation  for  large  slope  waves  (e.g.,  location  of  separation,  the  mean  flow  and  variance  in  the 
separation  zone).  It  is  amply  confirmed  that  the  LES  model  is  capable  of  predicting  these 
features  qualitatively  and  quantitatively.  The  results  from  this  study  have  been  published  in  the 
Journal  of  Fluid  Mechanics,  and  the  complete  paper  is  included  as  Appendix  A  of  this  report. 
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3.  TRACER  DISPERSION  RESULTS 


The  tracer  dispersion  calculations  require  the  solution  of  the  passive  scalar  conservation 
equations.  Since  the  source  size  will  often  be  less  than  the  grid  size  used  to  generate  the  velocity 
field,  Lagrangian  methods  are  used  for  the  scalar  field  so  that  subgrid  growth  can  be  computed. 
In  particular,  when  early-time,  near-source  behavior  is  of  primary  interest,  the  generalized 
Gaussian  puff  method  of  Sykes  and  Henn  (1995)  is  employed.  For  late-time  statistics,  i.e.,  when 
the  scalar  field  is  large  compared  to  the  length  scale  of  the  terrain,  a  particle  method,  which 
accounts  for  subgrid  diffusivity  with  a  random-walk  model,  is  used  since  only  the  overall 
centroid  and  second-moments  are  of  interest. 

This  study  focuses  on  dispersion  in  the  presence  of  simple,  idealized  terrain  so  that  the 
effects  of  parameters  such  as  slope  and  release  location  can  be  clearly  illustrated.  Given  the 
constraints  of  periodic  boundary  conditions  (on  the  velocity  field)  and  the  extensive  comparison 
with  velocity  measurements  over  sinusoidal  waves,  two-dimensional  sinusoidal  terrain  is  also 
used  for  the  dispersion  calculations.  Although  obviously  an  idealization,  this  is  a  reasonable 
prototype  of  ridge- valley  terrain. 

A  limited  range  of  terrain  slopes  and  atmospheric  stabilities  are  considered  in  the  two 
conference  papers  (Henn  and  Sykes  1996,1998).  Specifically,  the  maximum  slopes  are  0,  0.25 
and  0.5;  with  a  wavelength  of  2  km,  these  correspond  to  maximum  ridge  elevations  of  0,  159  m 
and  318  m,  respectively.  For  all  cases,  the  boundary  layer  depth  is  approximately  1000  m.  The 
maximum  elevations  are  a  substantial  fraction  of  this,  so  it  is  expected  that  the  terrain  will 
substantially  modify  the  boundary  layer  flow.  Two  uniform  surface  temperature  fluxes  are 
considered:  0.03  °Kms'1  and  0.0045  °Kms'1,  which  give  convective  velocity  scales  of  1  ms’1  and 
0.5  ms'1,  respectively.  The  higher  heat  flux  results  in  flow  which  is  dominated  by  large-scale 
convective  eddies,  while  the  lower  value  is  closer  to  neutral  stability  in  character.  The  LES 
domain  is  a  4  km  cube,  but  the  periodic  nature  of  the  velocity  field  is  utilized  to  allow  the 
dispersion  calculation  to  extend  beyond  the  nominal  velocity  field  domain.  Independent 
dispersion  realizations  are  obtained  by  varying  the  horizontal  location  of  the  source  within  the 
same  LES  velocity  field;  this  procedure  is  repeated  for  a  number  of  independent  velocity  fields. 
Thus,  ensemble  dispersion  statistics  are  generated  from  12  to  24  realizations. 
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The  principal  results  of  this  study  concern  the  sensitivity  of  short-range  dispersion  to 
release  location  as  well  as  the  enhanced  dispersion  over  longer  ranges.  Clearly  the  relatively 
large  terrain  slopes  considered  here  induce  complex  flow  features  such  as  re-circulation  zones, 
intermittent  separation,  vertical  mixing  and  vertical  shears,  which  increase  sensitivity  to  release 
location.  This  is  particularly  evident  with  the  lower  heat  flux  case;  in  contrast,  the  higher  heat 
flux  produces  stronger  vertical  mixing,  reducing  the  overall  effect  of  the  terrain.  Long-range 
calculations  also  clearly  point  out  the  effects  of  stronger  vertical  mixing.  For  instance,  the 
intermediate  terrain  slope  and  larger  heat  flux  results  in  no  enhanced  long-range,  streamwise 
dispersion  compared  with  the  flat  case,  whereas  the  low  heat  flux  case  gives  an  effective 
diffusion  coefficient  twice  the  flat  value.  Both  heat  fluxes  produce  significant  enhanced 
dispersion  with  the  largest  slope,  but,  as  expected,  the  low  heat  flux  shows  the  larger  increase 
(about  50%  greater  than  the  larger  heat  flux). 


6 


4.  CONCLUDING  REMARKS 


The  ability  of  LES  to  reproduce  the  observed  flow  features  in  boundary  layer  flow  over  a 
sinusoidally  varying  surface  has  been  demonstrated  through  comparison  with  laboratory  data. 
The  moderate  Reynolds  numbers  achieved  in  the  laboratory  are  sufficiently  high  for  the 
significant  turbulence  effects  to  be  similar  to  those  for  very  high  Reynolds  number  flow. 
Extensive  comparison  with  a  variety  of  experiments  and  observed  flow  statistics,  in  addition  to 
numerical  sensitivity  studies,  provides  sound  evidence  that  the  LES  flow  fields  accurately  model 
many  of  the  significant  features  of  flow  over  ridges.  One  of  the  most  striking  aspects  of  the  LES 
results  is  a  large  increase  in  lateral  velocity  fluctuation  intensity  on  the  upslope  of  the  ridges. 
The  increase  appears  to  be  associated  with  vortex  stretching  as  the  flow  is  accelerated  over  the 
crest  of  the  ridge,  but  the  lateral  velocity  component  was  not  measured  in  any  of  the  laboratory 
experiments,  so  this  feature  is  still  to  be  confirmed 

Further  LES  calculations  have  been  presented  for  buoyancy-driven  turbulence  over 
sinusoidal  ridges.  These  simulations  are  more  representative  of  atmospheric  flow,  where 
buoyancy  forcing  is  always  present.  Stability  categories  ranging  from  near-neutral  to  strongly 
convective  have  been  considered.  The  neutral  laboratory  comparisons  provide  confidence  that 
the  LES  can  properly  represent  the  flow  features,  since  buoyancy-driven  turbulence  is  more 
amenable  to  LES.  This  is  because  the  buoyancy-driven  eddies  are  large-scale  features,  and 
consequently  easier  to  resolve  in  the  numerical  calculation. 

Dispersion  calculations  for  a  passive  tracer  emitted  from  a  small,  instantaneous  source, 
located  at  several  positions  relative  to  the  terrain  variations,  were  then  presented.  The  dispersion 
results  clearly  show  the  increased  rate  of  diffusion,  particularly  in  the  streamwise  direction, 
associated  with  the  terrain-induced  turbulence  and  mean  flow  modifications.  The  results  also 
show  that  increasing  the  surface  buoyancy  flux  reduces  the  diffusion  rates,  due  to  the  increased 
vertical  mixing  rates.  The  major  mechanism  for  streamwise  extension  of  the  cloud  is  the 
increased  vertical  wind  shear  over  the  valley  regions;  increasing  the  vertical  mixing  reduces  the 
effective  shear.  An  extreme  example  of  the  shear  occurs  for  a  separated  flow  region,  where 
material  entrains  and  detrains  from  the  recirculation  zone  where  mean  the  streamlines  are  closed. 
Under  unstable  conditions,  the  exchange  rate  between  the  recirculation  zone  and  the  free  stream 
above  is  much  higher  than  under  neutral  conditions.  Thus  under  neutral  conditions,  we  see  much 
more  separation  between  particles  trapped  within  the  separation  region  and  particles  outside. 


7 


This  study  has  demonstrated  the  capability  of  LES  to  model  turbulent  dispersion  in 
boundary  layer  flow  over  terrain.  We  have  restricted  the  calculations  to  highly  idealized  terrain 
forms,  since  our  understanding  of  the  full  complexity  of  atmospheric  situations  can  only  be 
improved  in  stages.  LES  has  been  shown  to  reproduce  the  important  features  of  the  distorted 
turbulent  flow,  and  several  aspects  of  the  dispersion  process  have  been  elucidated,  using  the  LES 
velocity  fields.  Further  work  is  necessary  to  extend  the  calculations  to  other  terrain  forms. 
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Large-Eddy  Simulation  of  Flow  Over  Wavy  Surfaces 

Abstract 

Large-eddy  simulation  is  used  to  investigate  fully  developed  turbulent  flow  in  a  neutral 
channel  wherein  the  lower  wall  is  sinusoidal.  The  numerical  results  are  compared  with 
experimental  observations  for  wave  slopes  ranging  from  0  to  0.628.  Particular  emphasis  is 
placed  on  the  separated  flow  induced  by  a  large-amplitude  wave.  A  detailed  comparison  with 
the  data  of  Buckles  et  al.  (1984)  shows  generally  good  agreement.  LES  surface  pressures  are 
integrated  to  calculate  form  drag  as  a  function  of  wave  slope.  Drag  is  found  to  increase 
quadratically  with  slope  for  small-amplitude  waves,  with  a  somewhat  slower  increase  for  larger 
amplitudes.  However,  comparison  with  experimental  measurements  is  confounded  by 
uncertainties  with  the  values  reported  in  the  literature.  An  interesting  feature  characteristic  of  all 
wavy  surface  simulations  is  an  increase  in  transverse  velocity  fluctuations  on  the  wave  upslope. 
Although  the  precise  mechanism  responsible  is  not  known,  analysis  shows  it  to  be  associated 
with  temporally  persistent  vortex-like  structures  localized  near  the  surface.  The  magnitude  of  the 
fluctuation  increase  appears  to  scale  quadratically  with  slope  for  small-amplitude  waves,  in 
contrast  to  the  streamwise  fluctuations,  which  increase  linearly. 

1.  Introduction 

The  use  of  large-eddy  simulations  to  model  complex  flow  configurations  is  becoming 
increasingly  common.  One  of  the  earliest  applications  to  go  beyond  the  horizontally 
homogeneous  channel  case  was  in  modeling  turbulent  convective  flow  over  sinusoidal  terrain, 
e.g.  Krettenauer  and  Schumann  (1992),  Walko,  Cotton  and  Pielke  (1992)  and  Dombrack  and 
Schumann  (1993).  Simulations  over  flat  homogeneous  terrain  have  demonstrated  that  this  type 
of  flow  is  particularly  amenable  to  LES  since  the  dominant  thermally-driven  eddies  are  on  the 
scale  of  the  mixed-layer  height  (or  channel  depth)  and  so  are  typically  well-resolved  (Schmidt 
and  Schumann  1989;  Mason  1989).  While  LES  has  shown  that  wavy  terrain  can  modify  the 
large-scale  flow  in  the  convective  boundary  layer,  horizontally-averaged  profiles  are  not  very 
different  from  the  flat  case,  even  in  the  presence  of  a  mean  wind.  This  is  because  strong  vertical 
mixing  tends  to  reduce  the  velocity  shears  caused  by  even  large  surface  slopes.  However,  the 
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effects  of  terrain  on  more  neutrally  stable  flow  are  likely  to  be  much  more  dramatic  since  even 
relatively  small  slopes  can  induce  significant  variations  in  the  mean  flow  and  turbulence 
structure  (Sykes  1980;  Taylor  1977;  Belcher  and  Hunt  1998). 

This  research  effort  is  oriented  towards  atmospheric  flow  over  terrain  which,  while  idealized, 
is  nonetheless  realistic  in  terms  of  slope  and  wavelength.  However,  before  making  this 
extension,  it  is  desirable  to  verify  the  LES  against  laboratory  measurements.  Gong  et  al.  (1996) 
recently  presented  experimental  observations  and  LES  calculations  for  wind-tunnel  flow  over 
sinusoidal  terrain  with  a  growing  boundary  layer.  However,  the  flow  is  complicated  by  the  non¬ 
uniformity  of  the  boundary  layer  and  may  also  be  influenced  by  the  lateral  walls  of  the  wind 
tunnel.  Hanratty  and  coworkers  have  presented  an  extensive  set  of  measurements  of  channel  flow 
over  wavy  surfaces,  e.g.,  Zilker,  Cook  and  Hanratty  (1977),  Zilker  and  Hanratty  (1979),  Buckles, 
Hanratty  and  Adrian  (1984,  BHA  hereafter;  1986),  Frederick  and  Hanratty  (1988),  Kuzan, 
Hanratty  and  Adrian  (1989)  and  Hudson,  Dykhno  and  Hanratty  (1996).  We  will  focus  in 
particular  on  the  data  of  BHA,  who  present  extensive  measurements  of  the  mean  and  fluctuating 
velocity  for  turbulent  channel  flow  over  an  aerodynamically  smooth,  large-amplitude  sinusoidal 
surface.  However,  we  will  also  consider  smaller-amplitude  waves  in  order  to  establish 
confidence  in  the  model  as  well  as  examine  the  effects  of  varying  wave  slope.  Although  the  data 
we  are  comparing  with  are  all  for  smooth  walls  and  hence  Reynolds  number  dependent,  the 
Reynolds  numbers  are  high  enough  so  that  the  flows  are  fully  turbulent.  Nonetheless,  the  LES 
model  does  account  for  finite  Reynolds  number  effects  since  a  no-slip  velocity  boundary 
condition  is  used. 

Analytic  flow  studies,  e.g.,  Sykes  (1980),  Hunt  et  al.  (1988)  and  Belcher  et  al.  (1993),  have 
shown  that  terrain-induced  turbulence  modifications  are  subject  to  rapid  variation  close  to  the 
surface.  The  equilibrium  log-layer  relations  are  valid  only  very  near  the  boundary,  and  an  inner 
layer  develops  in  which  advection  and  distortion  of  the  turbulence  are  important,  with  a  vertical 
scale  much  smaller  than  the  horizontal  wavelength.  The  inner  layer  is  typically  about  5%  of  the 
wavelength,  making  explicit  resolution  of  eddies  in  this  region  very  difficult.  These  difficulties 
also  affect  the  ability  to  represent  surface  forces  correctly,  since  the  forces  are  closely  coupled  to 
the  details  of  the  flow  near  the  boundary.  However,  the  total  drag  force  is  also  related  to  the  total 
energy  dissipation  in  the  flow,  and  we  expect  wavy  surfaces,  especially  in  the  case  of  large-scale 
separation,  to  generate  resolvable  eddies  well  suited  to  modeling  by  LES. 

Direct  numerical  simulations  of  channel  flow  over  wavy  surfaces,  e.g.,  MaaB  and  Schumann 
(1994,1996)  ,  De  Angelis  et  al.  (1997)  and  Cherukat  et  al.  (1998)  have  successfully  reproduced 
many  of  the  laboratory  observations  for  moderate  Reynolds  number,  while  also  providing  great 
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detail  on  the  flow  structure.  For  instance,  both  De  Angelis  et  al.  and  Cherukat  et  al.  note 
increased  transverse  fluctuations  on  the  wave  upslope,  which  we  have  observed  for  a  range  of 
wave  amplitudes.  However,  the  application  of  DNS  is  limited  since  resolution  requirements 
become  prohibitive  for  high  Reynolds  number  flow. 

The  paper  is  organized  as  follows.  First  the  LES  model  is  described,  with  particular  emphasis 
on  features  relevant  to  finite  Reynolds  number  flow  over  wavy  surfaces,  since  the  infinite 
Reynolds  number,  Cartesian  version  of  the  model  has  been  described  elsewhere.  The  numerical 
methods  and  boundary  conditions  are  then  briefly  outlined.  After  summarizing  the  basic  flow 
configurations  and  parameters,  model  results  are  compared  with  a  number  of  experiments.  We 
first  present  results  from  flat  channel  calculations  and  then  proceed  to  examine  a  small-amplitude 
wavy  surface  case.  Next,  LES  results  are  compared  in  detail  with  the  data  of  BHA  for  separated 
flow  over  a  large-amplitude  wave.  A  limited  comparison  between  flows  over  waves  with 
different  amplitudes  is  then  made;  the  dependence  of  form  drag  on  wave  slope  is  examined  in 
particular.  Finally,  the  interesting  feature  of  increased  transverse  turbulent  fluctuations  on  the 
wave  upslope  is  investigated. 


2.  LES  model 


2.1  Model  equations 

The  basic  LES  model  has  previously  been  applied  to  buoyancy-driven  turbulence  (Sykes  and 
Henn,  1989;  Henn  and  Sykes  1992)  and  a  neutral  boundary  layer  (Sykes  and  Henn,  1992).  For 
the  current  application,  the  terrain-following  coordinate  transformation  of  Clark  (1977)  is  used 
for  accurate  representation  of  the  surface  boundary  conditions.  The  equations  of  motions  are 
based  on  the  filtered  Navier-Stokes  equations  and  follow  the  approach  of  Moin  and  Kim  (1982) 
in  explicitly  including  viscosity  so  that  a  no-slip  boundary  may  be  imposed.  The  resolved-scale 
Cartesian  velocity  components  =  (u,v,w)  are  solved  as  functions  of  transformed  coordinates 
x,.  =  (x,y,z),  defined  in  terms  of  Cartesian  coordinates  xi  =  (x,y,z)  by  x  =  x,y  =  ymd 
z  =  t](x,y,z).  The  transformation 

n=(z-h)/j  (i) 

where  J  =  1  -h/H ,  maps  the  domain  between  the  lower  wavy  surface  at  z  =  h  and  the  upper  flat 
channel  wall  at  z-H  onto  a  rectangular  domain.  The  lower  surface  is  a  two-dimensional  infinite 
wavetrain  defined  by  h(x,y )  =  a[l  +  cos(2^c//l)],  with  wavelength  A  and  amplitude  a. 
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Cartesian  derivatives  for  an  arbitrary  scalar  field  <p  are  given  in  terms  of  the  transformed 
coordinates  as 


dxj  d 
dxj  dxj 


0  -  G11  —r  <p 


dxj 


(2) 


where  the  summation  convention  is  assumed.  Note  that  J_1  =  det^G'-7  ).  Since  a  fluid  of  constant 
density  p  is  considered  here,  it  is  convenient  to  define  V  =  pj .  Then,  following  Krettenauer  and 
Schumann  (1992),  the  model  equations  are  given  concisely  as 


dx; 


(VGijUj)  =  0 


(3) 


Jr  (v“')+ar(VG;'"‘“') 3  'W^VGl‘p^^VGi‘T!‘^  vW{yG‘lD^+ VFA  (4) 


where p  is  the  dynamic  pressure,  Ttj  is  the  subgrid  stress  tensor,  vis  the  kinematic  viscosity  and 
Px  is  a  uniform  pressure  gradient  in  the  jc-direction  that  drives  the  flow.  It  should  be  noted  that 
this  is  not  intended  as  a  low  Reynolds  number  model;  the  viscous  term  is  only  important  in  a  thin 
layer  near  the  wall  and  is  included  principally  to  impose  the  no-slip  boundary  condition. 


The  subgrid  turbulence  model  uses  a  turbulent  kinetic  energy  transport  equation  based  on  the 
second-order  closure  model  of  Lewellen  (1977).  The  subgrid  velocity  variance,  q2,  is  obtained 
from  the  conservation  equation 


^(V)+ j~(va“ V) = -^—(vo’Fjyv 3 


3X; 


dX; 


a 


Gij^zr{VGkJq2) 
ox 


-2bV^~  (5) 
A 


where  is  the  subgrid  flux,  A  is  the  subgrid  turbulence  length  scale  (or  filter  scale)  and  the 
velocity  deformation  tensor  is  given  by 

a 


VD ..  =  • 

"  a; 


■(VG\  +  VGuUj) 


(6) 


The  first  term  on  the  right-hand  side  of  (5)  represents  subgrid  turbulence  production  due  to 
gradients  in  the  resolved  velocity  field,  the  next  two  terms  represents  transport  due  to  diffusive 
fluxes  and  the  last  term  models  the  dissipation.  The  subgrid  fluxes  are  given  by 

Vtt  =  S.qAVDt  (7) 

VF^SHqKj^(VG'q‘)  (8) 

The  various  empirical  model  constants  are  taken  from  Lewellen  (1977)  so  that  a  =  0.125, 
Sm=l/4andS„=l/3. 
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(9) 


Finally,  the  subgrid  length  scale  is  specified  algebraically  by 

1  _  i  1 

a2  ■  e +  a2 

where  l  is  a  function  of  distance  normal  to  the  nearer  of  the  top  or  lower  bounding  surfaces  and 
Amax  is  problem  dependent.  The  philosophy  of  Mason  and  Callen  (1986)  is  followed  in  setting 
the  subgrid  filter  scale  independent  of  the  numerical  grid,  subject  to  the  requirement  that  the 
resulting  large-scale  eddies  are  adequately  resolved.  While,  in  principle,  one  could  use  an 
arbitrarily  fine  grid  for  a  fixed  filter  scale,  in  practice  the  idea  is  to  reduce  the  length  scale  more 
or  less  in  proportion  with  the  grid  size  so  that  smaller-scale  turbulent  eddies  will  be  resolved. 
The  wall  function  t  is  simply  set  proportional  to  normal  distance  in  infinite-Reynolds  number 
applications.  However,  it  was  found  that  the  length  scale  must  be  damped  more  rapidly  as  the 
wall  is  approached  when  imposing  a  zero-slip  boundary  condition  for  finite-Reynolds  number 
flow.  Therefore,  following  Moin  and  Kim  (1982),  we  use  the  Van  Driest  exponential  damping 
function  such  that 

l  -  orn(l-exp(-n+/A+))  (10) 

where  n  is  distance  normal  to  the  nearest  wall,  a  =  0.65  from  Lewellen  (1977),  A+=25, 
n+  =  hTq/2/v  and  r0  is  the  local  instantaneous  wall  shear  stress. 

2.2  Numerical  implementation 

The  numerical  scheme  is  similar  to  that  used  in  the  Cartesian  model  (Sykes  and  Henn,  1989), 
although  the  coordinate  transformation  introduces  a  number  of  new  terms  and  other 
complications;  details  of  the  implementation  generally  follow  Clark  (1977).  The  model  uses  a 
standard  second-order  accurate,  finite  difference  energy-conserving  scheme  on  a  staggered  grid 
with  leapfrog  time  differencing.  A  small  amount  of  smoothing  (1%)  is  used  to  couple  the  two 
time  levels  and  prevent  the  time-splitting  instability.  Diffusion  terms  are  modeled  using  central 
differencing  with  the  DuFort-Frankel  approximation  to  provide  stable  integration.  The  DuFort- 
Frankel  scheme  is  second-order  accurate  and  consistent  with  the  advection-diffusion  equation  if 
d  =  K(kt/kz2)  is  small,  where  K  is  the  total  diffusivity.  At  the  wall,  K  =  v  and  d  is  less  than 
0.05  for  all  cases.  In  fact,  d  is  generally  less  than  0.01  over  most  of  the  domain  since  the 
timestep,  which  is  limited  by  the  vertical  Courant  condition  near  the  wall,  is  small.  The  grid 
spacing  is  uniform  horizontally  but  varies  vertically,  with  fine  resolution  at  the  walls  expanding 
to  a  section  .of  uniform  spacing  comparable  to  the  horizontal  spacing. 
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The  elliptic  equation  for  the  pressure  field  is  obtained  by  applying  the  finite-difference 
approximation  of  the  divergence  operator  (3)  to  the  partially  advanced  velocity  fields,  i.e., 
including  all  terms  from  (4)  except  the  pressure  gradients.  The  pressure  field  is  then  computed 
from  the  resulting  equation  such  that  zero  divergence  of  the  velocity  field  is  maintained. 
However,  the  non-orthogonal  coordinate  transformation  makes  the  equation  non-separable,  so 
the  direct  solver  of  the  Cartesian  model  must  be  used  in  an  iterative  scheme  as  in  Clark  (1977) 
and  Krettenauer  and  Schumann  (1992).  The  pressure  equation  is  written  in  the  form 

L0(p)  =  N{p)  +  D(ui)  (11) 

where  L0  is  the  finite-difference  representation  of  the  Cartesian  Laplace  operator  (i.e.,  the  case 
with  h  =  0),  N  represents  the  terms  of  the  operator  resulting  from  the  coordinate  transformation 
and  D  is  the  numerical  divergence  operator,  here  applied  to  the  partially  advanced  velocity  fields 
(denoted  by  a  tilde).  The  solution  is  obtained  using  the  direct  Cartesian  solver  for  each  iteration 
of  the  system 

A,(p("“>)  =  W(pw)  +  D(«,)  (12) 

where  the  superscripts  denote  iteration  and  the  previous  timestep  value  is  used  as  the  initial  guess 
p(0).  The  iteration  continues  until  the  maximum  change  in  p  is  less  than  10“5  times  the  maximum 
value  of  p  in  the  domain.  The  rate  of  convergence  depends  on  the  wave  slope  since  N  is 
generally  proportional  to  the  square  of  the  slope.  For  the  maximum  slope  of  0.628  considered 
here,  convergence  requires  5  to  7  iterations. 

2.3  Boundary  conditions 

Periodic  boundary  conditions  are  imposed  at  all  lateral  boundaries.  A  no-slip  condition  is 
imposed  on  the  lower  wall  by  setting  u  and  v  below  the  surface  equal  and  opposite  to  the  values 
on  the  first  level  above  the  surface.  For  some  of  the  flows  described  in  the  next  section,  a  no-slip 
condition  is  also  imposed  on  the  upper  wall.  However,  since  the  flow  details  at  the  top  of  the 
channel  are  generally  not  of  primary  interest,  the  upper  boundary  is  usually  treated  as  a  rough 
wall.  In  this  case,  a  logarithmic  profile  is  used  to  relate  the  top  wall  shear  stress  to  the  tangential 
velocity  half  a  grid  length  below  the  wall. 

The  momentum  fluxes  at  the  lower  surface  are  determined  by  the  no-slip  condition  at  the  wall 
and  by  assuming  flow  parallel  to  the  wall  at  the  first  grid  level.  The  total  momentum  flux  is 
given  by 


where  us  -  (u  ,v  ,hu+hy r)  is  the  velocity  vector  located  half  a  grid  length  above  the  lower 
surface,  hx  and  hy  are  the  surface  slopes  in  the  x-  and  y-di  recti  on  s,  respectively,  and  zn  is  the 
normal  distance  to  the  wall.  Note  that  for  the  purpose  of  setting  the  wall  stress,  w  is  defined  in 
(13)  so  that  u-  is  parallel  to  the  lower  surface.  The  normal  distance  is  related  to  the  transformed 
vertical  grid  spacing  at  the  wall,  Az0,  by 


= 


2 


(14) 


In  a  local  coordinate  system  defined  by  the  unit  surface  normal  hi  and  the  unit  vector  in  the 


-i 


direction  of  the  tangential  velocity,  =  uf  u ■  ,  the  full  stress  tensor  at  the  surface  is 


T,  =  Tr 


0  0  1 
0  0  0 
1  0  0 


(15) 


where  the  prime  denotes  the  local  coordinate  system.  The  stress  tensor  (15)  must  be  rotated  into 
the  global  Cartesian  frame  since  the  momentum  equations  use  the  Cartesian  velocity 
components.  Using  the  fact  that  hy  =  0,  this  rotation  results  in 


=  7A?0/2 


—2  h  u 

X  s 

—h  v 

X  5 

M. -h 


-hxvs  (l  -h])us 

0 

vr  2  h,.ur 


(16) 


which  defines  the  lower  boundary  condition  for  the  total  stress  tensor,  t ..  +  vD- . 


The  pressure  boundary  condition  is  set  to  insure  consistency  between  the  numerical 
discretizations  of  the  momentum  equation  (4)  and  the  pressure  source  term  N  from  (11)  when 
imposing  the  constraint  of  zero  normal  flow  at  the  wall.  By  requiring  that  both  ut  and  ui  (the 
partially  advance  velocity  used  in  (11))  satisfy  this  constraint,  we  obtain  a  numerical  expression 
of  the  form  Szp+  f(pz)  =  0  for  z  =  0,  where  8Z  denotes  a  vertical  finite  difference  and  / is  a 
function  of  the  vertically- averaged  pressure  pz.  The  numerical  expression  for  the  pressure 
source  term  is  given  by  N  -8zf  +NH(8.pz'),  where  Nh  is  a  horizontal  difference  operator. 
Thus,  setting  8zp  =  0  and  /  =  0  at  z  =  0  provides  all  the  necessary  boundary  conditions  for  the 
pressure  solver  while  maintaining  consistency  with  the  zero  normal  flow  condition. 


3.  Flow  parameters 

The  principal  focus  of  this  investigation  is  the  configuration  of  BHA,  i.e.,  channel  flow  over 
an  infinite  array  of  two-dimensional  large-amplitude  sinusoidal  waves.  However,  since  this  LES 
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model  had  not  previously  been  applied  to  flows  at  finite-Reynolds  numbers,  a  flat  channel  with 
Re  ~  12000  was  first  simulated  to  test  model  performance.  Then,  to  test  the  model  with  the 
transformed  equations,  channel  flow  over  a  small-amplitude  wavy  surface  was  calculated. 
Additionally,  intermediate  wave  amplitude  configurations  were  examined  as  part  of  a  limited 
investigation  into  variations  of  surface  pressure,  drag,  turbulence  intensity  and  flow  separation 
with  wave  slope. 


3.1  Numerical  experiments 

The  basic  flow  geometry  is  illustrated  in  Figure  1  and  all  configurations  are  summarized  in 
Table  1.  For  the  wavy  surface  cases,  the  standard  computational  domain  contains  two 
wavelengths;  test  calculations  with  a  single  wave  (not  shown)  gave  quite  different  results 
whereas  calculations  with  four  waves  were  very  similar  to  the  standard  domain  results.  This  can 
be  seen  from  the  values  for  case  BHA1L  given  in  Tables  3  (form  drag)  and  Table  4  (variance 
perturbations)  which  agree  within  10%  of  the  standard  domain  case  BHA1.  Note  also  that  grid 
resolutions  for  the  2a/ X  =  0,  0.031  and  0.2  cases  are  varied  along  with  the  subgrid  length  scale 
to  demonstrate  numerical  convergence. 

The  vertical  resolution  requirements  are  fairly  stringent  near  the  wall,  requiring  that  the 
lowest  grid  level  be  in  the  viscous  sublayer  below  z+  =  10  (using  local  wall  scaling).  Using  the 
experimental  results  to  estimate  the  maximum  mean  wall  shear  stress,  it  was  determined  that  grid 
levels  below  z/H  -0.005  and  0.01  were  required  for  the  BHA  and  flat  configurations, 
respectively.  The  minimum  vertical  grid  spacings  shown  in  Table  1  are  less  than  these  estimates 
to  accommodate  fluctuations  in  local  wall  stresses  resulting  from  very  thin  viscous  shear  layers. 

All  cases  were  initialized  from  previous  fully-turbulent  runs  and  continued  until  an 
approximate  equilibrium  between  the  imposed  constant  streamwise  pressure  gradient  and  the 
wall  stresses,  including  pressure  forces,  was  attained.  The  statistics  for  all  cases  presented  below 
are  from  periods  after  this  equilibrium  was  reached.  Averaging  was  done  over  a  time  of 
30  H/Ub  or  longer,  as  well  as  the  lateral  coordinate  and  position  of  equal  phase  relative  to  the 
surface  wave.  (The  flat  cases  were  averaged  over  the  lateral  and  streamwise  coordinates.)  The 
velocity  scale  used  here  is  given  by  a  bulk  velocity  that  is  independent  of  streamwise  location: 

Ub=(H-af'ffidz  (17) 

where  the  overbar  denotes  averaging.  Henceforth,  it  will  be  assumed  that  all  velocities  are 
normalized  with  respect  to  Ub- 
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Case 

L 

W 

X 

2  a/X 

Grid  Size 

A?m,n 

^^max 

^max 

Re 

F 

2 

1 

oo 

0 

48x48x90 

0.002 

0.02 

0.005 

11810 

FR 

2 

1 

oo 

0 

80x80x90 

0.002 

0.02 

0.003 

11820 

SI 

2 

1 

1 

0.031 

48x48x91 

0.002 

0.02 

0.005 

6560 

SIR 

2 

1 

1 

0.031 

80x80x102 

0.002 

0.02 

0.003 

5720 

S2 

1.984 

0.992 

0.992 

0.050 

48x48x71 

0.00124 

0.0248 

0.00620 

11050 

S3 

2.034 

1.017 

1.017 

0.100 

48x48x71 

0.00127 

0.0254 

0.00636 

10580 

S4 

2 

1 

1 

0.125 

48x48x67 

0.0025 

0.025 

0.00625 

10800 

BHA1 

2.143 

1.072 

1.072 

0.200 

48x48x71 

0.00134 

0.0268 

0.00670 

10600 

BHA1R 

2.143 

1.072 

1.072 

0.200 

80x80x99 

0.00134 

0.0167 

0.00402 

10450 

BHA1L 

4.286 

1.072 

1.072 

0.200 

96x48x71 

0.00134 

0.0268 

0.00670 

12000 

BHA2 

2.143 

1.072 

1.072 

0.200 

48x48x71 

0.00134 

0.0268 

0.00670 

5990 

BHA3 

2.143 

1.072 

1.072 

0.200 

48x48x71 

0.00134 

0.0268 

0.00670 

20060 

Table  1.  Non-dimensional  parameters  for  the  various  LES  calculations,  listed  in  order  of 
increasing  slope.  L,  W,  and  X  are  the  streamwise  length  (x-direction),  transverse  width  (y - 
direction)  and  wavelength,  respectively;  Azmm  and  Azmax  are  the  minimum  (at  the  lower  wall) 
and  maximum  vertical  transformed  grid  sizes,  respectively;  Amax  is  the  maximum  subgrid  filter 
length  scale;  and  the  Reynolds  number  is  defined  by  Re  =  UbS/v,  where  S  is  the  channel  mean 
half-height  and  Ub  is  the  bulk  velocity,  defined  in  the  text.  All  lengths  are  normalized  by  the 
mean  channel  depth,  H  -  a.  The  computational  grid  sizes  are  given  as  the  number  of  grid  points 
in  the  streamwise,  transverse  and  vertical  directions,  respectively. 


To  check  that  the  flow  was  sufficiently  stationary  for  computing  steady,  long-time  averages, 
the  time  evolution  of  Ub  and  the  total  pressure  and  shear  stress  forces  at  the  walls  were 
monitored.  As  an  illustration,  Figure  2  shows  the  time  history  of  the  integrated  surface  pressure 
force,  defined  in  equation  (19),  for  case  BHA1.  High-frequency  (period  ~H/Ub )  and  low- 
frequency  (period  ~  20 HjUh)  fluctuations  are  apparent.  Nonetheless,  it  is  evident  that  a 
meaningful  long-time  average  can  be  defined  over  the  averaging  period  shown  on  the  figure. 

All  simulations  were  performed  on  a  Windows  NT  workstation  employing  a  DEC  Alpha 
21164  300  MHz  processor.  For  the  standard  BHA1  case  (48x48x71  grid),  a  single  step  takes 
about  7  cpu-s.  Thus,  with  a  non-dimensional  timestep  of  8.6xl0~4,  the  integration  to  a  time  of  75 
took  approximately  170  cpu-h.  By  comparison,  the  flat  channel  case  F  (48x48x90  grid)  requires 
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about  5.5  cpu-s  per  step;  with  a  timestep  of  3.7x10  3  and  an  integration  time  of  166,  the  total  run 
time  was  about  69  cpu-h. 

It  will  be  noted  that  the  Reynolds  numbers  in  Table  1  do  not  match  the  experimental  values 
exactly.  Since  Px  is  held  constant  throughout  the  integration,  the  final  bulk  velocity  and,  hence, 
Reynolds  number  cannot  be  controlled  directly,  but  results  from  the  balance  between  wall 
stresses  and  imposed  pressure  gradient. 

In  presenting  the  LES  results,  the  resolved  turbulence  fluctuations  are  defined  as  u-UJ  where 
u"=  ut  -u.  while  the  total  fluctuations  include  the  subgrid  components,  u'u'j  =  u'u"-  .  We  also 

define  the  r.m.s.  streamwise  velocity  fluctuations  as  <7U  -  V«/2 ,  and  similarly  for  v  and  w.  Note 

that  in  presenting  results  for  the  flat  channel  and  small-amplitude  wave,  wall  scaling  is  denoted 
by  a  superscript  i.e.,  velocities  are  normalized  by  u*  =  ^(r0)  and  heights  by  v/w* ,  where 

angle  brackets  indicate  averaging  over  a  wavelength. 

3.2  Laboratory  experiments 

As  mentioned  in  the  introduction,  Hanratty  and  coworkers  at  the  University  of  Illinois  have 
presented  an  extensive  set  of  measurements  of  channel  flow  over  wavy  surfaces;  those  referred 
to  in  this  paper  are  summarized  in  Table  2.  The  experiments  cover  a  wide  range  of  slopes  and 
Reynolds  numbers  since  the  nature  of  the  flow  depends  strongly  on  these  parameters,  as 
illustrated  in  the  flow  regime  map  presented  by  Kuzan  et  al.  (1989).  We  have  concentrated  on 
experiments  with  Re  ~  104  where  the  mean  flow  separates  for  2ajX  somewhat  greater  than  0.05. 

All  experiments  were  conducted  in  a  rectangular  water  channel  with  a  cross  section  of  5.08  x 
61  cm.  A  long  flat  section  preceded  a  sinusoidal  wavetrain  machined  into  the  lower  channel  wall 
such  that  the  mean  cross  section  was  unchanged.  The  flat  section  was  long  enough  that  the  flow 
was  fully  turbulent  before  encountering  the  first  wave  and  the  channel  wide  enough  for  the  mean 
flow  to  be  essentially  two-dimensional.  Spatially  periodic  flow  was  established  fairly  rapidly; 
measurements  were  typically  taken  over  the  eighth  wave  in  a  series  of  ten.  (The  measurements 
of  Hudson,  et  al.  (1996)  were  taken  over  the  31st  of  36  waves.) 
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Reference 

Measurements 

Re 

2a/A 

Observed  Features 

Frederick  and 

Mean  and  fluctuating 

6400 

0.031 

Linear  shear  stress  response 

Hanratty  (1988) 

streamwise  velocity 
(LDV) 

Hudson,  et  al 

Mean  and  fluctuating 

3380 

0.100 

Separated  flow 

(1996) 

streamwise  and  vertical 

velocities,  Reynolds 
stress  (LDV) 

Zilker  and 

Mean  streamwise 

29300 

0.05 

Instantaneous  reverse  flow. 

Hanratty  (1979) 

velocity,  wall  pressure, 

non-linear  stress  response 

wall  shear  stress  (split- 

14730 

0.125 

Separated  flow 

film  thermal  probes. 

30000 

0.125 

pressure  taps,  electro¬ 
chemical  surface  probes) 

30530 

0.200 

Buckles,  et  al. 

Mean  and  fluctuating 

12300* 

0.125 

Separated  flow 

(1984,1986) 

streamwise  velocity,  wall 

10700* 

0.200 

pressure  (LDV  and 
pressure  taps) 


Table  2.  Summary  of  relevant  experiments  of  flow  over  wavy  surfaces  conducted  in  the 
laboratory  of  Hanratty  and  coworkers.  The  measurement  techniques  are  given  in  parentheses. 
*These  Reynolds  numbers  have  be  modified  from  their  published  values  to  conform  with  the 
bulk  velocity  definition  (17).  See  Section  6  for  further  discussion. 

4.  Flat  channel 

We  first  present  results  for  a  flat  channel  at  Re  ~  12000  with  no-slip  boundary  conditions  at 
both  walls,  similar  to  the  LES  study  of  Moin  and  Kim  (1982).  Figure  3  shows  the  mean 
streamwise  velocity  profiles  for  cases  F  and  FR  compared  with  the  data  of  Hussain  and  Reynolds 
(1970).  The  calculated  profiles  are  in  reasonable  agreement  with  the  data,  with  only  a  slight 
underprediction  evident  in  the  outer  part  of  the  wall  layer  (z+  =  10  to  30).  The  profiles  show 
very  little  sensitivity  to  resolution,  an  indication  of  numerical  convergence.  LES  profiles  from 
only  the  lower  wall  are  plotted,  but  those  from  the  upper  wall  are  virtually  indistinguishable,  an 
indication  of  adequate  statistical  averaging. 
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Figure  4  shows  profiles  of  the  turbulence  intensities,  resolved  and  total,  in  the  vicinity  of  the 
wall  while  Figure  5  shows  the  profiles  over  the  lower  half  of  the  channel.  The  peak  values  near 
the  wall  match  the  experiments  well,  as  do  the  values  in  the  middle  of  the  channel.  The 
streamwise  intensity  is  somewhat  overpredicted  for  z+  ~  20  to  80,  since  the  LES  shows  a  more 
gradual  fall  off  from  the  near-wall  peak.  It  is  evident  that  the  vertical  intensity  prediction  is 
determined  principally  by  the  subgrid  model  near  the  wall;  r33  increase  too  rapidly  from  the  wall 
but  nonetheless  shows  good  agreement  overall.  It  can  be  seen  that  the  resolved  components 
increase  as  the  subgrid  filter  scale  is  reduced  but  the  total  intensities  are  not  significantly 
changed,  indicating  an  appropriate  transfer  of  energy  from  the  subgrid  to  the  resolved  component 
as  resolution  is  increased.  The  agreement  with  the  experimental  data  is  comparable  to  the  LES 
of  Moin  and  Kim  (1982). 

5.  Small-amplitude  wavy  surface 

Before  presenting  results  for  the  large-amplitude  wave,  we  simulate  flow  over  a  small- 
amplitude  wave  that  tests  the  coordinate  transformation  while  still  yielding  results  close  to  the 
flat  channel.  This  flow  is  of  interest  in  its  own  right,  so  we  will  compare  with  experimental  and 
theoretical  results.  This  will  enable  us  to  examine  the  model  predictions  without  the 
complicating  effects  of  separation. 


5.1.  Mean  velocity 

Results  from  cases  SI  and  SIR,  with  2ajX- 0.031,  are  compared  with  the  measurements  of 
Frederick  and  Hanratty  (1988).  Figure  6  shows  the  wavelength-averaged  mean  streamwise 
velocity.  The  LES  results  are  averaged  along  lines  of  constant  rj  while  the  data  are  averaged 
along  lines  of  constant  height  above  the  wave  surface.  For  this  small  slope,  the  difference 
between  the  two  averaging  definitions  is  negligible.  The  predicted  profiles  are  generally  in  close 
agreement  with  the  measurements  and  the  two  numerical  resolutions  result  in  nearly  identical 
profiles.  As  with  the  comparison  shown  in  Figure  3  for  the  flat  channel,  the  LES  profiles  are 
somewhat  low  over  the  approximate  range  10  <  z+  <  30,  although  the  discrepancy  is  slightly  less 
here.  The  flat  channel  data  of  Hussain  and  Reynolds  (1970)  are  also  shown  in  the  figure;  it  is 
evident  that  the  wavelength  average  matches  these  data  quite  closely,  as  would  be  expected  for 
linear  perturbation  at  this  small  slope. 

The  velocity  profiles  at  different  locations  along  the  wave  are  shown  in  Figure  7.  Profiles  are 
paired  such  that  they  should  be  symmetric  around  the  wavelength  average  (shown  as  a  dashed 
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line),  assuming  perfectly  linear  behavior.  The  overall  agreement  with  the  data  is  good, 
particularly  above  the  viscous  sublayer,  z+>10.  The  largest  discrepancies  occur  near  the 
surface  of  the  upslope  profiles  at  x/A  =  0.6  and  0.7,  where  the  measurements  show  greater  flow 
acceleration.  It  seems  that  the  LES  predictions  in  these  cases  show  more  linear  behavior,  i.e.,  the 
paired  LES  profiles  are  more  nearly  symmetric  about  the  wavelength  average  than  the 
experimental  data. 


In  a  study  of  flow  over  hills  using  second-order  closure  and  mixing-length  models,  Newley 
(1985)  examines  the  streamwise  velocity  perturbation  profile  above  the  wave  crest  and  finds  that 
the  maximum  can  be  expressed  as 

laX) 


A  w„ 


akul(£m) 


(18) 


where  u0  is  the  unperturbed  mean  streamwise  velocity,  ak  is  the  maximum  slope,  wavenumber 
k  =  2nj  A,  £m  is  the  height  of  the  maximum  perturbation  above  the  crest,  f0is  an  outer  layer  scale 
set  to  k'1  for  periodic  hills  and  J3  is  an  0(1)  constant.  A  fit  to  the  second-order  closure  results 
gives  ft  ~  1  (and  ft  ~  0.8  for  the  mixing-length  model).  The  LES  results,  using  the  wavelength 
average  velocity  for  u0,  give  /?=  1.5.  This  discrepancy  partly  reflects  the  fact  that  lm  is  higher 
with  the  LES  model:  £m/A  =  0.013  compared  with  approximately  0.006  for  Newley’s  second- 
order  closure  calculations  for  hills  with  similar  slopes.  Of  course,  Newley’s  calculations  are  for 
rough  surfaces,  so  these  results  are  not  directly  comparable,  especially  since  £m  is  just  barely 
above  the  viscous  sublayer.  If  we  re-examine  the  data  of  Figure  7  in  terms  of  the  perturbation 
velocity  near  the  surface,  as  is  done  in  Frederick  (1986),  we  see  that  the  prediction  of  £m  matches 
the  experiment  quite  well.  As  can  be  seen  from  Figure  8,  the  LES  profiles  at  the  wave  crest  and 
trough  show  good  overall  agreement  with  the  data,  although  the  overprediction  at  x/A  =  0.5, 
z+  >  10  is  more  noticeable  with  the  linear  scale.  The  maximum  perturbations  above  the  crest  are 
only  about  10%  higher  than  the  experimental  values,  indicating  that  the  LES  prediction  of  ft  is 
reasonably  consistent  with  the  experimental  data. 


5.2.  Turbulence  intensity 


Streamwise  turbulence  intensity  measurements  from  Frederick  (1986)  are  shown  in  Figure  9 
along  with  the  LES  results.  The  predictions  compare  favorably  overall  with  the  experiment  in 
that  the  basic  variations  with  location  are  reproduced.  The  turbulence  intensities  generally  are 
out  of  phase  with  the  pressure  gradients:  a  positive  (adverse)  pressure  gradient  has  a  destabilizing 
effect  and  hence  increases  the  velocity  fluctuations  while  the  opposite  is  true  for  negative 
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(favorable)  pressure  gradients.  Thus  the  peak  values  increase  from  profiles  at  x/A  =  0.1  to  0.5 
and  decrease  thereafter  to  x/A  =  1 .  The  LES  peak  values  match  the  measurements  quite  well 
regarding  both  magnitude  and  location  (typically  around  z+  ~  20).  Near  the  wall,  the  LES 
predictions  are  generally  low  and  show  greater  differences  between  the  paired  profiles.  Away 
from  the  wall,  z+  >  100,  the  LES  also  tends  to  underpredict  the  fluctuation  intensity.  However, 
flat  channel  measurements  from  the  same  facility  (Frederick,  1986)  show  10-20%  greater 
centerline  intensities  than  typically  given  in  the  literature,  e.g.,  Hussain  and  Reynolds  (1975). 

To  relate  the  LES  results  to  some  of  the  theoretical  models  of  flow  over  small-amplitude 
waves,  e.g.,  Sykes  (1980)  and  Newley  (1986),  we  show  the  streamwise  and  vertical  velocity 
variance  perturbation  fields  for  the  high  resolution  case  SIR  in  Figure  10.  (The  transverse 
component  will  be  discussed  in  Section  9.)  The  streamwise  component  is  defined  as 
Am'2  =  u'2-(u'2y  with  the  vertical  component  defined  similarly.  There  is  qualitative  agreement 

with  the  closure  model  predictions:  a  distinct  surface  layer  is  evident;  the  perturbations  maximize 
above  and  are  almost  180°  out  of  phase  with  the  surface  layer;  the  maximum  positive  streamwise 
perturbation  occurs  over  the  trough  while  the  largest  negative  perturbation  occurs  over  the  crest. 
A  quantitative  comparison  with  the  models  is  not  meaningful  due  to  differences  in  the  bottom 
boundary  conditions;  the  inner  layers  defined  in  Sykes  (1980)  and  Belcher  et  al.  (1993)  both 
would  fall  within  the  viscous  sublayer  of  the  LES,  thus  invalidating  their  assumption  of  a  log 
profile.  It  is  interesting  to  note  that  Am'2  above  the  surface  layer  is  out  of  phase  with  Aw'2 . 
Whereas  the  streamwise  variance  decreases  over  the  wave  crest,  the  vertical  variance  increases 
there,  a  result  not  predicted  by  closure  models  or  rapid  distortion  theory  but  observed  in 
measurements  over  hills  (Belcher  et  al.  1993).  It  should  also  be  noted  that  the  maximum  in  Am'2 
is  localized  close  to  the  surface  compared  with  Aw'2 .  Thus,  although  the  maximum  magnitude 
of  Am'2  is  nearly  ten  times  greater  than  Aw'2 ,  it  decays  rapidly  with  height  and  so  for  the  region 
r]/A  >  0.1,  the  perturbations  are  of  comparable  magnitude. 

The  results  shown  in  Figure  10  are  determined  to  some  extent  by  the  subgrid  model, 
especially  near  the  surface.  Therefore,  it  is  important  to  examine  how  much  of  the  total 
variances  are  explicitly  resolved  in  the  LES.  The  wavelength-averaged  fractions  of  resolved 
variances,  defined  as  ^u"2^ j (u'2^j  for  the  streamwise  variance  and  similarly  for  the  transverse 

and  vertical  variances,  are  shown  in  Figure  11  for  cases  SI  and  SIR.  It  can  be  seen  that  more 
than  90%  of  the  streamwise  variance  is  explicitly  resolved  over  the  entire  domain  for  both 
numerical  resolutions.  In  contrast,  only  10%  of  the  vertical  variance  is  resolved  very  close  to  the 
surface  in  SI  while  the  minimum  resolved  fraction  of  the  transverse  variance  approaches  60%. 
Higher  resolution  results  in  a  noticeable  improvement:  the  minimum  resolved  fraction  increases 
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to  35%  for  the  vertical  variance  and  75%  for  the  transverse  component.  In  the  standard 
resolution  case  SI,  all  components  are  more  than  80%  resolved  for  ri/A>  0.1.  This  level  of 
resolution  is  achieved  for  T]/ A  >  0.05  with  the  higher-resolution  of  SIR.  Thus,  the  variance 
perturbations  shown  in  Figure  10  are  essentially  resolved  above  the  surface  layer;  the 
calculations  within  the  surface  layer  must  be  viewed  with  some  caution  however. 

5.3.  Surface  shear  stress 

Zilker  (1976)  presents  surface  wall  stress  data  for  this  flow;  a  comparison  with  the  LES 
results  is  shown  in  Figure  12.  The  model  predictions  are  close  to  the  measurements  in  both 
phase  and  magnitude,  although  the  experimental  scatter  suggests  caution  in  making  any 
definitive  conclusions.  It  is  evident  that  the  wall  stress  data  show  a  phase  lead  relative  to  the 
wavecrest.  Fitting  the  data  with  a  two-harmonic  Fourier  series,  Frederick  (1986)  finds  a  phase  of 
51°  for  the  dominant  first  harmonic.  The  LES  shows  a  slightly  smaller  phase  shift  of 
approximately  45°.  These  results  are  consistent  with  theoretical  models,  e.g.  Hunt  et  al.  (1988) 
and  Sykes  (1980)  and  reflect  the  fact  that  the  maximum  perturbation  velocity  near  the  surface  is 
shifted  upstream  (Newley  1986).  The  shear  stress  perturbation  is  not  exactly  symmetric:  the 
increase  just  ahead  of  the  wave  crest  is  of  larger  magnitude  than  the  decrease  that  occurs  just 
upwind  of  the  trough.  Accordingly,  the  low  stress  region  around  the  trough  is  of  somewhat 
greater  extent  than  the  high  stress  region. 

6.  Separated  flow  over  a  large-amplitude  wavy  surface 

We  now  turn  to  the  BHA  configuration,  2 a/ A  =  0.2,  that  exhibited  a  large  separated  region 
extending  over  most  of  the  wave  trough.  There  is  some  uncertainty  in  the  definition  of  the  bulk 
velocity  used  in  BHA.  They  give  a  definition  that  is  a  function  of  position,  which  seems  like  an 
unnecessary  complication.  On  the  other  hand.  Buckles  (1983)  uses  a  bulk  velocity  with  a 
nomalizing  depth  equal  to  the  minimum  channel  depth,  H  -  2a,  rather  than  the  mean  channel 
depth.  That  is  also  consistent  with  integrals  of  the  tabulated  velocity  profiles  given  in  Buckles 
(1983)  and  used  here  for  comparison  with  the  LES.  Therefore,  for  this  section  only,  we  use  a 
definition  of  bulk  velocity  given  by 

Ub=(H-2a)~lj”udz 

Note  that  using  this  definition  gives  a  Reynolds  number  12%  higher  compared  with  that  obtained 
with  the  bulk  velocity  definition  (17).  This  discrepancy  is  reflected  in  the  Reynolds  numbers 
given  in  Tables  2  and  3. 
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6.1.  Mean  velocity  and  streamlines 


Profiles  of  the  mean  streamwise  velocity  for  the  simulations  BHA1  and  BHA1R  are  compared 
with  the  experimental  data  in  Figure  13.  Clearly,  both  the  measurements  and  the  LES  show 
strongly  separated  flow.  The  predicted  magnitude  of  the  reversed  flow  in  the  recirculation  zone, 
roughly  between  xj  A  =  0.1  and  0.7,  tends  to  be  slightly  greater  than  observed,  but  overall  the 
profile  variations  along  the  wave  are  faithfully  reproduced.  Some  locations  show  an 
underprediction  in  the  middle  of  the  channel,  but  this  appears  to  be  due  to  variations  in  the 
measurements.  We  note  that  the  profiles  reported  by  BHA  contain  discrepancies  of  up  to  3%  in 
the  integrated  mass  flux,  which  is  necessarily  conserved  in  the  LES. 

Profiles  at  the  crest,  trough  and  two  intermediate  positions  are  shown  on  a  log  scale  in  Figure 
14  to  emphasize  the  near  surface  region.  The  LES  profiles  (only  BHA1R  is  shown  since  BHA1 
is  very  similar)  are  in  reasonable  agreement  with  the  measured  data.  In  particular,  they  match 
well  over  the  crest  except  for  a  slight  underprediction  very  close  to  the  surface.  Significant  shear 
can  be  seen  very  close  to  the  surface  at  the  crest.  Above  this,  the  velocity  profile  is  nearly 
constant  for  a  short  while  before  increasing  to  the  channel  centerline  maximum.  The  LES 
overpredicts  the  reverse  flow  in  the  trough,  but  shows  better  agreement  on  the  up-  and 
downslope  locations.  Note  that  the  LES  profiles  show  slightly  negative  velocity  at  x/A  =  0.7 
while  the  measurements  are  slightly  positive,  indicative  of  a  somewhat  further  downstream 
reattachment  point  in  the  LES. 

The  mean  streamwise  velocity  field  and  streamlines  are  shown  in  Figure  15  for  BHA1R. 
(The  results  for  BHA1  are  similar.)  The  strong  shear  over  the  crests  is  clearly  shown,  but  it 
should  also  be  noted  that  there  is  still  significant  shear  centered  more  or  less  around  the  VF  =  0 
streamline.  The  recirculation  zone,  bounded  by  the  T  =  0  streamline  and  the  wave  surface,  is 
somewhat  deeper  than  that  shown  in  BHA  (their  Figure  5);  the  bounding  streamline  has  a 
maximum  elevation  approximately  90%  of  the  wave  height  compared  with  60%  in  BHA. 
Although  the  LES  velocity  profiles  agree  well  with  the  observations  over  the  bulk  of  the 
recirculation  zone,  the  overprediction  of  the  reverse  flow  near  the  surface  results  in  the  integral 
defining  the  streamlines,  namely,  ¥  =  \u  dz' ,  going  to  zero  at  a  greater  height. 

A  detailed  examination  of  the  LES  and  experimental  data  shows  some  discrepancy  in  the 
precise  extent  of  the  recirculation  region.  Separation  (in  the  mean  flow)  occurs  at  x/A  =  0.14  in 
the  experiment  compared  with  0.06  for  both  BHA1  and  BHA1R;  reattachment  occurs  at  0.69  in 
the  experiment  and  at  0.73  and  0.75  for  BHA1  and  BHA1R  respectively.  Since  the 
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experimenters  note  that  velocity  measurements  are  inaccurate  very  close  to  the  surface  (below 
z/ / l  =  2.6xl0“3),  there  is  some  uncertainty  in  the  precise  locations  of  the  mean  separation  and 
reattachment  points.  Nonetheless,  this  discrepancy  has  little  adverse  effect  on  predicting  the 
salient  features  of  the  mean  flow. 

6.2.  Intermittency  and  flow  variability 

It  is  suggested  in  BHA  that  the  location  of  the  zero  mean  streamwise  velocity  contour 
corresponds  closely  with  the  loci  of  50%  instantaneous  reverse  flow,  i.e.,  y=  0.5  where  the 
intermittency  yis  defined  as  the  fraction  of  time  that  u  is  positive.  This  is  clearly  confirmed  for 
the  LES  in  Figure  15,  where  the  y=  0.5  line  originates  very  close  to  the  separation  point  and 
terminates  very  close  to  the  point  of  reattachment,  as  in  BHA.  The  y=  0. 1  contour  indicates  that 
there  is  a  significant  fraction  of  reverse  flow  at  the  top  of  the  mean  recirculation  zone,  i.e.,  at 
approximately  the  height  of  the  wave  crests.  The  location  of  this  line  matches  well  that  shown  in 
BHA,  although  they  note  that  it  falls  surprisingly  high  above  the  recirculation  zone.  This  is  not 
the  case  with  the  LES  due  to  its  deeper  recirculation  zone. 

The  intermittency  contours  imply  large  variability  in  the  flow  associated  with  the  recirculation 
zone.  Plots  of  instantaneous  velocity  fields,  examples  of  which  are  shown  in  Figure  16,  indicate 
that  large-scale  flow  features  are  responsible  for  this  variability.  Clearly,  flow  separation  results 
in  large-scale  “flapping”  of  the  shear  layer  as  well  as  other  complex  behavior.  For  example,  the 
character  of  the  flow  is  strikingly  different  in  each  of  the  four  troughs  shown  in  Figure  16.  In 
one  case,  reverse  flow  fills  most  of  the  “valley”  above  the  trough  while,  in  another,  reverse  flow 
is  confined  to  a  shallow  region  near  the  wave  surface.  In  a  third  case,  reverse  flow  seems  to  have 
been  ejected  away  from  the  surface.  Finally,  we  see  a  number  of  small  areas,  typically  on  the 
downslope,  where  the  reverse  flow  appears  to  have  separated  so  that  there  is  forward  flow  very 
near  the  surface.  It  is  clear  from  these  examples  that  the  mean  streamlines  are  not  at  all 
representative  of  the  instantaneous  flow,  but  they  do  illustrate  the  region  where  recirculating 
flow  is  likely  to  be  observed. 


6.3.  Higher-order  statistics 

The  relative  streamwise  turbulent  fluctuations  are  shown  in  Figure  17,  where  it  can  be  seen 
that  the  LES  results  are  in  close  agreement  with  the  data.  Only  the  total  fluctuation  intensities 
are  plotted  since  the  subgrid  components  are  almost  negligible  everywhere.  BHA1R  profiles  are 
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in  close  agreement  with  those  of  the  lower  resolution  BHA1;  there  are  slight  differences  in  some 
details,  but  statistics  for  the  two  resolutions  differ  by  less  than  a  few  percent. 

The  highest  fluctuation  intensities  are  associated  with  the  free  shear  layer  that  detaches  from 
the  surface  at  the  separation  point  and  is  elevated  above  the  separation  zone.  The  magnitude  of 
this  maximum  intensity  remains  fairly  constant  over  much  of  the  wave.  A  reduction  is  seen  after 
the  flow  reattaches  and  speeds  up  as  it  approaches  the  wave  crest.  These  features  can  be  seen 
more  clearly  in  Figure  18,  where  turbulence  intensity  fields  are  shown  for  the  three  coordinate 
directions.  The  locus  of  maximum  streamwise  intensity  extends  over  most  of  the  trough  at  a 
height  equal  roughly  to  the  maximum  wave  height,  before  elevating  somewhat  towards  the  crest. 
This  locus  approximately  coincides  with  the  midpoint  of  the  vertical  shear  region  in  the  mean 
velocity  field  (see  Figure  15(a)).  The  vertical  extent  of  this  free  shear  layer  is  evident  from  the 
“elbows”  in  the  mean  velocity  profiles  of  Figure  13,  particularly  for  x/A  =  0.3  and  0.5.  The 
magnitude  of  the  maximum  streamwise  intensity  increases  rapidly  from  a  value  of  about  0.20 
very  close  to  the  crest  to  around  0.25  a  quarter  of  the  wavelength  downstream. 

As  discussed  in  BHA,  remnants  of  the  free  shear  layer  separating  from  the  upstream  wave  are 
manifested  as  an  elevated  (z/H  ~  0.25  to  0.35)  region  of  increased  streamwise  turbulence 
intensity  over  the  wave  crest.  It  can  be  viewed  simply  as  a  continuation  of  the  locus  of 
maximum  intensity  mentioned  above  though  a  periodic  domain.  However,  this  elevated 
maximum  decays  rapidly  downstream  of  the  crest. 

A  moderate  increase  in  the  vertical  fluctuation  intensity,  crw,  above  the  wave  crest  may  also  be 
associated  with  the  shear  layer  remnant.  However,  the  maximum  intensity  occurs  over  the 
trough  and  is  probably  related  to  the  large-scale  eddies  responsible  for  the  variations  of  the 
separated  flow,  as  illustrated  in  Figure  16. 

Unlike  au  and  aw,  the  transverse  fluctuation  intensity,  <rv,  does  not  show  any  particular  feature 
associated  with  the  free  shear  layer.  However,  there  is  a  very  pronounced  maximum  near  the 
surface  on  the  wave  upslope  whose  effects  extend  over  the  downstream  crest.  This  is  a  feature 
which  seems  characteristic  of  all  wavy  surface  cases  that  we  have  examined.  Since  it  is  rather 
unexpected  and,  to  our  knowledge,  has  not  been  observed  experimentally,  we  will  discuss  it 
further  in  Section  9. 

The  contour  plots  of  streamwise  and  vertical  turbulence  intensities  in  Figure  18  appear 
qualitatively  similar  to  those  shown  in  Hudson  et  al.  (1996)  for  separated  flow  over  a  smaller 
wave  (la/ A  =  0.1)  at  a  lower  Reynolds  number  (Re  =  3380).  In  fact,  the  locations  and  even  the 
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magnitudes  of  the  maximum  observed  intensities  are  reasonably  close  to  the  LES  predictions. 
Hudson  et  al.  correlate  the  different  locations  of  <JU  and  crvv  maxima  with  the  different  direct 
production  terms  (involving  the  Reynolds  stresses  and  mean  velocity  gradients).  They  also  point 
out  that  the  lag  in  vertical  intensity  maximum  relative  to  the  streamwise  intensity  is  consistent 
with  the  initial  direct  production  of  streamwise  fluctuations  (through  u'w'du/dz  and  u'2  du/dx) 
and  the  subsequent  transfer  of  energy  into  vertical  fluctuations  via  the  pressure-strain  term. 

Hudson  et  al.  also  show  contours  of  the  Reynolds  stress  -u'w' ,  which  again  is  in  qualitative 
agreement  with  the  LES  results  shown  in  Figure  18.  It  can  be  seen  that  the  elevation  of  the 
maximum  Reynolds  stress  is  very  close  to  that  of  maximum  <ru  and  is  clearly  associated  with  the 
separated  shear  layer.  The  negative  values  close  to  the  surface  on  the  upslope  are  an  artifact  of 
the  Cartesian  coordinate  system  (Hudson  1993);  they  assume  positive  values  if  rotated  into  a 
boundary-layer  or  streamline  coordinate  system. 

Profiles  of  skewness  and  flatness  of  the  streamwise  velocity,  defined  as  S  =  n"3/cr3  and 
F  =  u"A  /  <rAu,  respectively,  are  compared  with  the  BHA  data  for  x/ X  =  0.3  in  Figure  19.  The 
agreement  between  the  LES  profiles  and  the  measurements  is  encouraging  considering  the 
difficulty  in  obtaining  reliable  higher-order  statistics.  Away  from  the  wall,  the  profiles  do  not 
differ  greatly  from  those  computed  or  measured  previously  over  a  flat  wall,  e.g.,  Moin  and  Kim 
(1982).  The  mid-channel  values  of  S  ~  -0.4  and  F  ~  3  are  close  to  the  flat  channel  values  and 
result  from  velocity  fluctuations  produced  by  slow-moving  fluid  arriving  from  the  wall  region. 
However,  the  profiles  clearly  differ  from  flat  channel  flow  near  the  wall.  For  instance,  S 
decreases  from  its  maximum,  which  is  associated  with  the  elevated  shear  layer,  to  a  small 
negative  value  at  the  wall.  In  contrast,  the  flat  channel  profile  of  skewness  maximizes  at  the 
wall.  This  is  also  true  for  flatness.  However,  for  the  BHA  flow,  F  shows  some  rather  complex 
behavior  close  to  the  wall,  particularly  around  (z  -  h)/H  ~  0.05.  While  not  inconsistent  with  the 
measurements,  it  remains  to  be  seen  if  this  is  real  or  an  artifact  of  the  LES  model.  The  LES 
profiles  do  show  some  sensitivity  to  resolution,  especially  for  F  near  the  wall.  This  is  not 
surprising  since  only  the  resolved  velocity  fluctuations  are  considered. 

Contour  plots  of  skewness  and  flatness  are  shown  in  Figure  20.  It  can  be  seen  that  a  region  of 
high  skewness  occurs  just  downstream  of  the  crest,  reaching  its  maximum  around  the  separation 
point.  The  locus  of  the  maximum  skewness  seems  to  follow  approximately  along  the  middle  of 
the  recirculation  zone.  Large  negative  values  of  skewness  are  seen  very  close  to  the  wave 
surface  over  the  region  of  reverse  flow  while  near-surface  positive  skewness  seems  to  be 
associated  with  attached,  forward-moving  fluid.  Both  observations  are  consistent  with  the  view 
that  velocity  fluctuations  near  the  wall  are  caused  by  fast-moving  fluid  being  mixed  down  from 
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above.  But,  the  skewness  is  then  necessarily  negative  in  the  reverse  flow  region.  The  region  of 
high  flatness  is  closely  correlated  with  the  region  of  high  skewness,  suggesting  that  the  process 
of  mixing  faster-moving  fluid  to  the  wall  is  intermittent. 

6.4.  Surface  pressure  and  shear  stress 

The  final  comparisons  with  BHA  are  of  the  pressure  and  shear  stress  on  the  wavy  surface.  An 
accurate  prediction  of  these  quantities  is  clearly  an  important  aspect  of  LES  since  they  determine 
the  integrated  forces  on  the  wave,  which  is  of  interest  in  aeronautics,  oceanography,  meteorology 
and  other  fields.  Variations  of  the  pressure  force  with  slope  are  discussed  in  Section  8;  here  we 
focus  on  the  BHA  configuration. 

Figure  21  shows  distributions  of  non-dimensional  wall  pressure  pw  and  shear  stress  Tw 
compared  with  the  experimentally  measured  distributions,  where 

P»  =  (y)  and  r"  =(y) 

The  pressure  plots  have  been  adjusted  so  that  the  LES  and  observations  match  at  the  wave 
trough.  It  can  be  seen  that  both  BHA1  and  BHA1R  predict  the  pressure  increase  on  the  upslope 
quite  well.  However,  the  magnitude  of  the  pressure  trough  at  the  wave  crest  is  under-predicted 
in  both  cases  and  some  sensitivity  to  resolution  is  evident.  The  wall  shear  stress  is  generally  well 
predicted  although  the  phasing  of  the  peak  stress  lags  somewhat  compared  to  the  experiment. 
Note  that  the  experimental  surface  shear  stress  was  estimated  by  assuming  a  linear  velocity 
profile  from  the  lowest  measurement  point  to  the  wall  and  there  is  some  uncertainty  associated 
with  these  measurements. 

The  total  forces  on  the  lower  wall  per  unit  area  in  the  x-direction  can  be  calculated  by 
integrating  the  curves  in  Figure  21.  Defining  the  non-dimensional  pressure  and  shear  stress 
forces  (per  unit  length)  as 


the  experimental  values  are  Fx  =  0.0261  and  Tx  =  0.0028.  This  compares  with  Fx  =  0.0253  and 
0.0265  for  BHA1  and  BHA1R,  respectively;  the  corresponding  results  for  the  shear  stress  are 
Tx  =  0.0016  and  0.0014.  The  accurate  prediction  of  the  pressure  force,  commonly  called  the 
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form  drag,  is  due  almost  entirely  to  the  fact  that  the  pressure  rise  in  Figure  21  matches  the  data 
well.  The  under-prediction  of  the  pressure  trough  is  not  critical  to  the  integrated  force  in  the  x- 
direction  since  the  local  pressure  force  there  acts  vertically.  However,  the  predicted  shear  stress 
force  is  significantly  lower  than  the  experimental  observations  because  the  biggest  discrepancy 
occurs  on  the  upslope  (which  is  opposite  to  what  we  just  noted  for  the  pressure  distribution). 
This  discrepancy  stems  from  an  underprediction  of  the  acceleration  on  the  upslope  and  the 
resulting  delayed  flow  reattachment.  Although  the  difficulty  in  measuring  velocity  close  to  the 
surface  may  be  a  factor  in  the  discrepancy,  in  any  case,  the  magnitude  of  the  shear  stress 
component  is  small  compared  with  the  pressure  contribution,  so  the  total  drag  force  is  reasonably 
well  predicted  by  the  numerical  simulation. 

7.  Variations  with  slope 

The  LES  model  has  been  used  to  examine  the  effects  of  varying  wave  slope.  Figure  22  shows 
mean  streamwise  velocity  profiles  for  2a/ A  =  0.05,  0.10  and  0.20  (cases  S2,  S3  and  BHA1  in 
Table  1).  It  is  evident  that  the  mean  flow  separates  for  2a/ A  =  0.1  but  the  0.05  case  is  marginal 
in  this  regard.  Other  profile  locations  (not  shown)  between  x/A  =  0.3  and  0.5  reveal  that  case  S2 
does  show  a  very  small  region  of  flow  separation.  This  finding  is  consistent  with  the  flow 
regime  map  of  Kuzan  et  al.  (1989),  which  also  indicates  that  the  higher  slope  cases  will  show 
clearly  separated  mean  flow.  The  depth  of  the  shear  region  above  the  trough  clearly  scales  with 
the  wave  amplitude.  This  is  also  true  for  the  height  of  maximum  streamwise  turbulence 
intensity,  shown  in  Figure  23.  The  turbulence  intensity  increases  almost  linearly  with  slope, 
especially  away  from  the  wave  surface. 

Figure  24  shows  the  variation  in  surface  pressure  and  shear  stress  with  wave  slope.  The 
pressure  curves  show  a  definite  change  in  character  as  slope  increases.  The  small-amplitude 
wave  produces  a  smooth,  slightly  asymmetric,  curve  that  is  similar  to  the  one  shown  in  Zilker 
and  Hanratty  (1979)  for  this  slope,  but  with  a  somewhat  smaller  peak.  However,  there  may  be 
some  Reynolds  number  effect  given  that  the  measurements  by  Zilker  and  Hanratty  are  for 
Re  =  30000.  As  slope  increases,  a  nearly  constant  pressure  region  over  the  separated  flow 
becomes  more  pronounced.  The  shape  of  the  2a/ A  =  0.2  curve  is  similar  to  that  in  Zilker  and 
Hanratty  (as  well  as  BHA),  but  apparently  the  LES  significantly  underpredicts  the  total  pressure 
perturbation  (peak-to-trough).  (As  we  saw  above,  the  2a/ A  =  0.2  case  shows  good  overall 
agreement  with  the  data  of  BHA  if  the  predictions  and  observations  are  matched  in  the  region  of 
uniform  pressure  at  the  trough.)  A  comparison  of  the  2 a/ A  =  0.1  case  with  the  2a/ A  =0.125  data 
in  Zilker  and  Hanratty  also  reveals  an  underprediction  of  the  pressure  perturbation.  In  this  case, 
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though,  the  curves  are  also  qualitatively  different.  As  noted  above,  the  LES  curve  shows 
evidence  of  separation  by  a  relatively  flat  section.  However,  although  Zilker  and  Hanratty 
observed  separation  at  this  slope,  their  pressure  data  do  not  show  a  similar  feature. 

The  surface  shear  stress  curves  shown  in  Figure  24(b)  clearly  illustrate  the  regions  of 
separated  or  stagnant  flow  for  the  three  slopes.  The  shear  stress  for  2a/  A  =  0.05  is  nearly 
symmetric  (implying  very  little  contribution  to  the  increase  in  total  surface  force)  and  becomes 
slightly  negative  approximately  from  x/A  =  0.2  to  0.5.  In  light  of  the  velocity  profiles  in  Figure 
22,  it  is  evident  that  the  flow  here  is  only  weakly  separated  and  is  probably  better  characterized 
as  stagnant.  The  region  of  negative  shear  stress  clearly  increases  with  slope,  as  does  the 
magnitude  of  the  negative  stress.  The  local  maxima  both  2afA  =  0.1  and  0.2  cases  exhibit  around 
x/A  =  0.3  is  rather  surprising,  but  results  from  a  small  reduction  in  the  magnitude  of  the  reverse 
flow  there.  Since  measurements  of  surface  shear  stress  are  somewhat  unreliable,  we  cannot 
ascertain  if  this  feature  is  realistic.  Because  the  reattachment  point  (the  second  zero  stress 
crossing)  moves  downstream  with  increasing  slope,  the  subsequent  rate  of  increase  in  shear 
stress  tends  to  be  greater,  although  the  location  of  the  maximum  also  tends  to  move  downstream 
slightly. 

The  pressure  curves  in  Figure  24(a)  certainly  show  significant  variations  with  slope. 
However,  probably  of  greater  interest  than  the  curves  themselves  is  the  resulting  integrated  force 
on  the  wave  surface,  especially  as  a  function  of  slope.  This  is  examined  in  the  next  section. 


8.  Form  drag  calculations 

In  this  section,  we  present  LES  calculations  of  form  drag  for  a  number  of  wave  slopes;  see 
Table  1  for  a  summary  of  cases.  These  calculations  include  some  variations  with  resolution  and 
Reynolds  number  since  the  results  may  be  somewhat  controversial  and  we  want  to  establish  the 
sensitivity  and  consistency  of  the  numerical  calculations. 

Figure  25  shows  the  variations  in  form  drag  as  a  function  of  slope.  In  comparing  with  the 
data  of  Zilker  and  Hanratty,  as  is  done  in  Belcher  et  al.  (1993)  and  Gong  et  al.  (1996),  we  found 
that  the  LES  form  drag  predictions  were  approximately  twice  as  large  as  the  reported  values. 
However,  we  checked  the  Zilker  and  Hanratty  results  by  integrating  the  corresponding  pressure 
curves  as  presented  in  their  paper.  As  seen  from  Table  3,  the  numerically  integrated  forces  are 
around  twice  the  values  stated  in  Zilker  and  Hanratty  and  are  reasonably  consistent  with  the  LES 
results  as  well  as  the  measurements  of  BHA.  These  calculations  are  shown  in  Figure  25  as  the 
solid  symbols.  We  have  employed  this  method  on  the  data  of  BHA  and  get  very  close  agreement 
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with  the  reported  form  drag  (within  6%).  We  also  show  the  form  drag  derived  from  the  pressure 
data  presented  in  Buckles  et  al.  (1986)  for  la/ A  =  0.125  {Re  =  12300),  although  they  do  not  give 
their  own  calculated  value.  The  form  drag  in  this  case  appears  to  be  quite  high,  being  roughly 
equal  with  the  la/ A  =  0.2  value,  although  not  that  much  higher  than  the  numerically  integrated 
Zilker  and  Hanratty  result  for  Re  -  14730.  This  suggests  there  may  be  some  dramatic  change  in 
the  nature  of  the  flow  at  this  slope,  possibly  with  some  Reynolds  number  dependence,  or  it  may 
point  out  some  difficulties  and  uncertainties  in  the  measurements.  Note  that  the  measurements 
of  Gong  et  al.  are  still  in  reasonable  agreement  with  the  adjusted  values  of  Zilker  &  Hanratty, 
particularly  if  the  normalization  in  Gong  et  al  is  corrected  for  the  fact  that  the  centerline  velocity 
is  approximately  20%  larger  than  the  bulk  velocity.  However,  the  predictions  in  Belcher  et  al. 
appear  to  be  much  lower  than  the  adjusted  Zilker  and  Hanratty  forces  as  well  as  those  of  BHA 
and  the  LES  in  this  study. 

We  have  examined  sensitivity  to  resolution,  domain  length  and  Reynolds  number  for 
2 a/ A  =  0.2.  As  can  be  seen  from  the  figure  and  Table  3,  the  form  drag  calculated  from  the  LES 
appears  to  be  relatively  insensitive  in  the  range  of  Reynolds  numbers  considered.  This  is  in 
contrast  to  Zilker  and  Hanratty,  especially  for  the  numerically  integrated  pressure  forces, 
although  there  is  certainly  some  question  as  to  the  reliability  of  these  force  results. 

It  should  be  pointed  out  that  the  LES  results  are  subject  to  some  uncertainty  since  they  depend 
entirely  on  modeling  the  flow  near  the  surface  where,  as  has  been  shown,  the  explicit  resolution 
of  the  turbulent  motions  is  effectively  reduced.  However,  the  effective  resolution  increases  with 
slope  since  larger-scale  motions  are  generated  and  the  form  drag  predictions  are  self-consistent 
in  that  a  quadratic  dependence  on  slope  is  shown  over  a  range  of  amplitudes,  i.e.,  a  good  fit  is 
given  by 

Fx  =  0.12  (ak)2 

This  relationship  appears  to  be  valid  up  to  around  2a/ A  =  0.1;  it  starts  to  fall  off  as  slope  and 
hence  separation  increases.  This  behavior  is  consistent  with  theoretical  expectations,  e.g.  Wood 
and  Mason  (1993). 
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Reference 

2a/JL 

Re 

Fx 

^  1 
^  | 

Jo rw  dx 

Zilker  &  Hanratty  (1979) 

0.05 

29300 

1.26xl0-3 

2.89xl0"3 

0.125 

14730 

1.155xl0-2 

2.29x1 0-2 

0.125 

30000 

1.003xl0"2 

1.64xl0~2 

0.20 

30530 

1.37xl0~2 

2.56xl0~2 

Buckles  et  al.  (1986) 

0.125 

12300 

- 

3.08xl0"2 

BHA 

0.20 

10700 

3.27xl0“2 

3.46xl0~2 

LES -SI 

0.031 

6560 

- 

1.24xl0-3 

SIR 

0.031 

5720 

- 

1.31xl0~3 

S2 

0.05 

11050 

- 

2.90xl0“3 

S3 

0.10 

10580 

“ 

1.29xl0~2 

S4 

0.125 

10800 

- 

1.70xl0"2 

BHA1 

0.20 

10600 

- 

3.16xl0-2 

BHA1R 

0.20 

10450 

- 

3.34xl0"2 

BHA1L 

0.20 

12000 

- 

3.21xl0-2 

BHA2 

0.20 

5990 

- 

3.20xl0'2 

BHA3 

0.20 

20060 

3.19xl0“2 

Table  3.  Form  drag  for  various  wave  slopes  and  Reynolds  numbers.  Fx  is  the  normalized  form 
drag  given  in  the  cited  references;  the  last  column  is  the  form  drag  computed  using  published 
surface  pressure  data  (calculated  by  the  present  authors).  Note  that  the  values  for  BHA  and 
Buckles  et  al.  (1986)  are  scaled  by  (H  -  a)2  / (H  -  2a)2  so  they  are  consistent  with  £4  defined  by 
(17). 


9.  Increased  transverse  fluctuations 

As  shown  in  Figure  18,  the  transverse  fluctuations,  v'2 ,  show  a  marked  increase  on  the 
upslope  close  to  the  wave  surface.  In  fact,  the  largest  velocity  fluctuations  are  found  in  the 
lateral  component,  exceeding  even  the  streamwise  component  maximum,  located  in  the 
separating  shear  layer.  A  similar  feature  has  recently  been  reported  by  De  Angelis  et  al.  (1997) 
and  Cherukat  et  al.  (1998)  in  their  DNS  studies  of  flow  over  waves  with  smaller  slopes,  although 
not  as  pronounced  as  in  Figure  18.  The  magnitude  and  limited  spatial  extent  of  the  v'2  increase 
strongly  suggests  a  localized  energy  production  mechanism  associated  with  the  wave  slope  and 
we,  therefore,  investigate  the  lateral  velocity  fluctuations  in  more  detail. 
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An  examination  of  instantaneous  flow  fields  at  r/jX  =  0.025,  such  as  that  shown  in  Figure  26, 
indicates  structures  elongated  in  the  streamwise  direction;  they  are  most  easily  seen  in  the 
transverse  velocity  field.  However,  they  do  not  resemble  the  streaks  seen  in  flat  boundary  layers: 
these  are  broader  laterally  and  limited  in  streamwise  extent.  These  broad  areas  of  coherent 
transverse  velocity  are  also  evident  in  a  plot  of  the  secondary  flow  vectors  at  x/X  =  0.75,  Figure 
26(c).  There  are  intense  vortex-like  transverse  motions  close  to  the  surface,  but  they  do  not  seem 
to  be  clearly  defined  closed  circulations  or  streamwise  rolls. 

The  velocity  fluctuations  associated  with  these  intense  transverse  velocity  regions  are  strongly 
affected  by  the  presence  of  the  solid  boundary,  since  they  occur  so  close  to  the  wall.  This  most 
likely  is  the  reason  for  the  preferential  manifestation  of  the  kinetic  energy  in  the  v-component, 
rather  than  equally  amongst  v  and  vv.  For  a  streamwise  vortex,  an  idealized  reflective  boundary 
condition  at  the  wall  reduces  the  normal  component  of  the  velocity  to  zero  and  doubles  the 
transverse  velocity  component,  giving  a  strong  enhancement  of  v'2  at  the  surface.  The  wall 
layer  drag  reduces  the  velocity  very  close  to  the  surface,  but  the  large  increase  in  v'2  on  the 
lower  side  of  the  vortex-like  structures  is  clearly  evident  in  the  LES  results. 

The  persistence  of  these  structures  is  illustrated  in  Figure  27,  which  shows  the  transverse 
velocity  component  at  x/X  =  0.75  and  r)j/ 1  =  0.0125  as  a  function  of  time  and  transverse 
position.  Examples  from  BF1A1  and  the  flat  channel  case  F  are  shown.  It  can  be  seen  that  the 
flat  channel  contours  are  rather  isotropic  in  appearance,  indicating  an  intermittent  flow  structure 
that  is  not  strongly  correlated  in  time  or  along  the  transverse  direction.  In  contrast,  the  contours 
for  the  wavy  channel  midway  on  the  upslope  show  broad  streaks  of  positive  or  negative  velocity. 
Although  these  features  meander  laterally  and  the  magnitude  of  the  v-velocity  in  the  streak  is 
modulated  over  time,  the  identity  of  the  structure  is  maintained  throughout,  providing  strong 
evidence  of  their  stability. 

Both  De  Angelis  et  al.  (1997)  and  Cherukat  et  al.  (1998)  suggest  that  the  transverse  velocity 
fluctuation  increase  is  produced  by  the  large-scale  structures  associated  with  flow  separation  and 
reattachment.  This  is  not  necessarily  the  case.  Although  the  v'2  maximum  is  localized  near  the 
reattachment  point  over  the  large-amplitude  wave,  LES  of  flow  over  smaller  wave  slopes  also 
show  a  similar  feature  on  the  wave  upslope.  Figure  28  shows  the  transverse  fluctuation 
intensities  for  the  three  waves  considered  in  Section  8.  Note  that  the  maximum  on  the  upslope  is 
present  for  all  waves,  including  the  smallest  slope  that  exhibits  very  weak  separation.  (The  mean 
reattachment  point  in  that  case  is  at  x/X  =  0.5 ,  well  upstream  of  the  v'2  maximum.). 
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Figure  29  shows  the  transverse  variance  perturbation  for  case  SIR,  as  was  done  for  the 
streamwise  and  vertical  components  in  Figure  10.  It  can  be  seen  that  for  a  small-amplitude 
wave,  a  noticeable  increase  in  Av/2  on  the  upslope  is  evident,  but  it  is  accompanied  by  an  almost 
equal  decrease  on  the  downslope.  (However,  an  examination  of  the  larger-amplitude  waves 
reveals  that  this  near  anti-symmetry  breaks  down  as  slope  increase  so  that  the  magnitude  of  the 
upslope  (positive)  perturbation  becomes  much  larger  than  the  magnitude  of  the  downslope 
(negative)  perturbation.)  In  contrast  to  Figure  10,  there  is  little  qualitative  agreement  with 
closure  model  predictions. 

The  perturbation  maximum  in  Figure  29  is  close  to  the  surface  and  thus  is  influenced  by  the 
subgrid  model,  but  Figure  11  shows  that  at  least  75%  of  the  transverse  variance  is  resolved  at 
that  location.  Also,  the  fraction  of  resolved  variance  increases  with  slope  (not  shown),  so  results 
should  be  more  reliable  for  the  the  larger-amplitude  cases. 

The  results  shown  in  Figure  28  indicate  that  v/2  increases  with  wave  steepness.  To  quantify 
this,  we  examine  the  normalized  maximum  perturbations  for  the  streamwise  and  transverse 
variances  as  functions  of  wave  slope,  given  in  Table  4.  In  can  be  seen  that  A u'2  scales  closely 
with  slope  for  small-amplitude  waves,  in  accord  with  linear  theories  (Sykes  1980).  However,  the 
Av'2  increases  suggest  a  quadratic  dependence  on  slope,  with  a  fall  off  at  larger  slopes  similar  to 
the  form  drag. 


LES  Case 

ak 

Aw'2  /ak 

A  v'2/{akf 

SI 

0.098 

7.70xl0"2 

2.41X10'1 

SIR 

0.098 

8. 43x1  O'2 

2.50X10’1 

S2 

0.157 

7.91xl0'2 

1.93X10'1 

S3 

0.314 

9.13xl0‘2 

2.64x10'' 

S4 

0.393 

8.18xl0'2 

1.86  xlO"' 

BHA1 

0.628 

8.67xl0"2 

1.69x10' 

BHA1R 

0.628 

6.47xl0"2 

1.34x10-' 

BHA1L 

0.628 

9.55xl0'2 

1.81x10"' 

Table  4.  Maximum  variance  perturbations  for  LES  cases  defined  in  Table  1.  The  streamwise 
perturbations  are  normalized  by  the  wave  slope  ak;  the  transverse  perturbations  are  normalized 
by  the  square  of  the  wave  slope. 
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Examination  of  the  conservation  equation  for  v'2  shows  that,  since  there  are  no  mean 
gradients  in  the  lateral  direction,  direct  production  terms  are  zero  and  the  transverse  velocity 
fluctuations  must  be  driven  by  the  pressure-strain  term,  p'dv'/dy,  as  pointed  out  by  De  Angelis 
et  al.  (1997).  It  is  well-known  that  pressure-strain  terms  redistribute  energy  from  one  component 
to  another  and  this  is  clearly  the  case  with  v'2 :  energy  from  the  streamwise  and/or  vertical 
components  are  extracted  through  the  correlation  between  pressure  perturbations  and  transverse 
velocity  gradients.  The  proximity  of  the  wall  is  certainly  a  factor,  as  mentioned  above,  but  the 
basic  mechanism  that  gives  rise  to  this  correlation  is  unknown.  It  may  be  due  to  amplification  of 
some  flow  instability,  but  it  is  not  related  to  the  large-scale  shear  instability  described  in  Phillips 
et  al.  (1996)  and  possibly  observed  in  the  experiments  of  Gong  et  al.  (1996),  since  that 
mechanism  has  a  streamwise  scale  greater  than  the  wavelength.  The  increased  v'2  we  observe  is 
much  more  limited  in  extent.  A  possible  mechanism  is  the  Taylor-Gortler  instability  which 
produces  streamwise  vortices  such  as  observed  in  boundary  layer  flow  over  a  concave  surface, 
e.g.,  Tani  (1962),  So  and  Mellor  (1975)  and  Hoffman  et.  al.  (1985).  The  trough  region  presents  a 
concave  surface,  and  we  might  therefore  expect  production  of  coherent  streamwise  vortices  in 
the  local  flow.  However,  the  magnitude  of  the  v'2  increase  observed  in  the  LES  calculations 
appears  to  be  significantly  larger  than  experimental  observations  in  curved  boundary  layer  flow, 
but  this  may  be  due  to  enhancement  by  the  shear  distortion  in  the  flow  over  the  waves. 
Experimental  flows  have  generally  involved  minimal  streamwise  acceleration  in  an  attempt  to 
isolate  the  curvature  mechanisms.  In  contrast,  the  vorticity  production  in  flow  over  waves  occurs 
in  a  region  of  strong  acceleration  and,  hence,  vortex  stretching.  Streamwise  vorticity  is  therefore 
enhanced  as  the  flow  accelerates  toward  the  wave  crest,  increasing  the  v-velocity  and  also 
stabilizing  any  existing  vortices. 

10.  Conclusions 

Large-eddy  simulations  of  flow  over  a  wavy  surface  have  been  presented  and  compared  with 
available  experimental  data  at  Reynolds  numbers  of  about  104.  The  viscous  effects  are  only 
important  close  to  the  smooth  wall,  and  the  LES  grid  size  normal  to  the  wall  is  chosen  to  resolve 
the  viscous  sublayer.  A  range  of  wave  amplitudes  has  been  studied,  from  small  slopes  with 
almost  linear  flow  perturbations  up  to  steep  waves  with  large-scale  flow  separation.  Particular 
attention  has  been  paid  to  the  large  amplitude  wave  experiment  of  BHA. 

In  spite  of  the  fact  that  surface  elevation  variations  are  known  to  induce  complex  distortion  of 
the  turbulence  fields  close  to  the  surface,  the  moderate  resolution  of  the  LES  appears  sufficient  to 
capture  many  aspects  of  the  experimental  flow.  This  is  partially  due  to  the  fact  that  the 
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streamwise  velocity  fluctuations,  often  the  only  component  measured  in  the  laboratory,  are 
relatively  well  resolved  in  the  LES.  Nevertheless,  the  response  of  the  mean  flow  and  the 
streamwise  velocity  variance  to  different  wave  amplitudes  is  accurately  predicted  even  close  to 
the  surface. 

Small  amplitude  waves  illustrate  the  nature  of  the  perturbations,  and  confirm  the  turbulence 
closure  model  predictions  of  an  "inner  layer"  structure  between  the  wall-layer  and  the  outer  flow 
on  the  scale  of  the  horizontal  wavelength.  The  LES  results  show  rapid  variation  in  the  velocity 
variances  within  a  layer  of  about  5%  of  the  wavelength;  this  is  generally  consistent  with  closure 
theory  predictions  but  quantitative  comparison  is  precluded  by  the  log-layer  assumption  of  the 
analytical  models. 

At  large  amplitude,  the  flow  is  dominated  by  the  separating  shear  layer  which  generates  large- 
scale  eddies,  allowing  more  of  the  turbulent  energy  to  be  resolved  by  the  LES.  The  predicted 
separation  region  is  more  extensive  than  in  the  experiment,  but  the  turbulence  and  mean  flow 
fields  are  only  sensitive  to  this  discrepancy  very  close  to  the  wall.  The  predicted  structure  of  the 
streamwise  velocity  fluctuations,  including  the  higher  order  statistics  of  skewness  and  flatness,  is 
consistent  with  the  experimental  data. 

A  comparison  of  predicted  form  drag  with  experimental  observations  has  pointed  out  some 
uncertainties  in  the  observations.  In  particular,  it  appears  that  the  drag  coefficients  reported  in 
Zilker  and  Hanratty  (1979)  are  inconsistent  with  their  surface  pressure  distributions.  However, 
the  form  drag  reported  in  BHA  (which  is  consistent  with  their  surface  pressure  data)  is  close  to 
the  LES  value.  The  LES  results  are  subject  to  uncertainty  since  resolution  is  effectively  reduced 
near  the  surface,  but  they  are  self-consistent  in  that  the  form  drag  shows  a  quadratic  dependence 
on  slope  for  small-amplitude  waves,  as  predicted  by  analytical  models.  A  slight  fall  off  with 
slope  is  predicted  for  larger  amplitudes,  again  consistent  with  expectations,  e.g.,  Wood  and 
Mason  (1993). 

One  of  the  most  striking  features  of  the  numerical  calculations  is  a  dramatic  increase  in  lateral 
velocity  variance  in  a  very  localized  region  close  to  the  surface  on  the  upslope  of  the  wave.  This 
is  a  general  feature  of  all  our  wavy  surface  simulations.  The  increase  apparently  scales  with  the 
square  of  the  slope,  in  contrast  to  the  other  two  velocity  components,  which  scale  linearly  with 
slope  for  small  amplitude  waves.  The  eddies  responsible  for  the  high  variance  are  persistent 
structures,  elongated  along  the  upslope,  with  relatively  large  lateral  scale.  The  lateral  motions 
are  driven  by  pressure  variations,  constrained  by  the  tangential  plane  of  the  wall,  but  the  basic 
mechanism  for  their  generation  is  unclear.  This  phenomenon  has  been  observed  in  the  DNS 
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studies  of  De  Angelis  et  al.  (1997)  and  Cherukat  et  al.  (1998),  but,  as  far  as  we  are  aware,  not 
measured  experimentally,  and  further  investigation  is  warranted. 
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Figure  1.  Schematic  diagram  of  channel  flow  over  a  periodic  wavy  lower  wall.  The 
maximum  channel  depth  is  H;  the  LES  computational  domain  has  a  length  L.  For  most 
cases  with  a  wavy  surface,  L  is  twice  the  wavelength  A,  the  transverse  computational 
width,  W,  is  equal  to  A  and  A  ~  H .  For  the  flat  channel  cases,  L  =  2 H  and  W  -  H . 
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Figure  2.  Time-history  of  integrated  surface  pressure  force  or  form  drag,  defined  by  (19)  for  case 
BHA1.  The  averaging  period  used  to  generate  steady-state  statistics  is  shown  on  the 


Figure  3.  Mean  velocity  profile  for  flat  channel  cases  compared  with  experimental  data.  Solid 
lines  are  standard  resolution  case  F;  dashed  lines  are  higher  resolution  case  FR;  symbols 
denote  measurements 
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Figure  4.  Turbulence  intensities  for  cases  F  and  FR  compared  with  experimental  data.  The  solid 
lines  are  the  total  fluctuations;  dashed  lines  are  the  resolved  components.  The  high 
resolution  case  FR  is  shown  with  dark  lines,  (a)  Streamwise  turbulence  intensity  (b) 
transverse  turbulence  intensity,  (c)  vertical  turbulence  intensity. 
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Figure  5.  As  Figure  4,  but  height  above  wall  scaled  by  channel  half-depth,  (a)  Streamwise 
turbulence  intensity,  (b)  transverse  turbulence  intensity,  (c)  vertical  turbulence  intensity. 


Figure  6.  Mean  wavelength-averaged  velocity  profile  for  cases  SI  and  SIR  (2a/ A  =  0.031) 
compared  with  experimental  data.  The  standard  resolution  case  S 1  is  shown  as  the  solid 
line;  the  higher-resolution  case  SIR  is  shown  as  the  dashed  line;  symbols  denote 
experimental  data. 
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Figure  7.  Profiles  of  mean  streamwise  velocity  at  various  locations  along  the  surface  wave 
compared  with  data  of  Frederick  and  Hanratty  (1988).  For  the  downslope  locations 
x/X  =0.1  to  0.5,  thin  solid  lines  are  case  SI,  thin  dashed  lines  are  case  SIR  and  open 
circles  are  the  experimental  data.  For  the  upslope  locations  x/X  =  0.6  to  1.0,  thick  solid 
lines  are  case  SI,  thick  dashed  lines  are  case  SIR  and  solid  circles  are  the  experimental 
data.  The  wavelength  average  from  SIR  is  shown  with  the  long-dash,  short-dash  pattern. 
(The  curves  from  the  two  LES  cases  are  virtually  indistinguishable.) 
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Figure  8.  Profiles  of  mean  streamwise  velocity  perturbation  at  the  wave  crest  and  trough.  Solid 
lines  are  case  SI;  dashed  lines  are  SIR;  open  circles  are  the  experimental  data  of 
Frederick  (1986). 
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Figure  11.  Fractions  of  velocity  variances  resolved  by  the  LES.  The  standard  resolution  case  SI 
is  shown  with  solid  lines,  the  higher-resolution  case  SIR  with  dashed  lines.  The  lines 
labeled  “u”  are  fractions  for  the  streamwise  component  (defined  in  the  text);  “v”  and  “w” 
denote  the  transverse  and  vertical  components,  respectively. 
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Figure  12.  Non-dimensional  shear  stress  distributions.  Solid  line  is  case  SI;  dashed  line  is  SIR 
symbols  are  the  data  of  Zilker  (1976)  from  two  adjacent  waves. 
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Figure  13.  Profiles  of  mean  streamwise  velocity  at  various  locations  along  the  surface  wave 
compared  with  the  data  of  BHA.  Solid  lines  are  LES  case  BHA1;  dashed  lines  are 
BHA1R;  symbols  are  measurements. 
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Figure  14.  Profiles  of  mean  streamwise  velocity  using  a  logarithmic  vertical  axis.  Only  profiles 
from  BHA1R  are  shown  for  comparison  with  the  data.  (The  profiles  from  BHA1  are 
very  similar.)  The  lines  denote  LES  profiles,  symbols  denote  BHA  measurements. 
x/A  =  0.3:  short  dashes;  x.  x/A  =  0.5 :  long  dashes;  o.  x/A  =  0.7 :  long-dash,  short-dash; 
□  .  x/A  - 1:  solid  line;  +. 
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Figure  17.  Profiles  of  mean  streamwise  turbulence  intensity  at  various  locations  along  the  surface 
wave  compared  with  the  data  of  BHA.  Solid  lines  are  BHA1;  dashed  lines  are  BHA1R; 
symbols  are  measurements.  Total  LES  turbulence  intensities  are  shown. 


BHA1;  dashed  lines  are  BHA1R;  symbols  are  data  of  BHA 
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Figure  22.  Profiles  of  mean  streamwise  velocity  for  3  slopes:  case  BHA1  (2a//l  =  0.2),  solid 
lines;  case  S3  (2a//l  =  0.1),  long  dashes;  case  S2  (2a/ A  =  0.05),  short  dashes. 


Figure  23.  Profiles  of  mean  streamwise  turbulence  intensity  for  cases  BHA1,  S3  and  S2.  Lines 
patterns  as  in  Figure  22. 


Figure  25.  Form  drag  vs.  wave  slope  as  given  in  Table  3.  O:  data  of  Zilker  and  Hanratty  (1979); 
•  :  integration  of  Zilker  and  Hanratty  pressure  data;  □:  data  of  BHA;  ■  :  integration  of 
BHA  pressure  data;  A:  integration  of  LES  surface  pressure.  Dashed  line  is  proportional 
to  (akf. 
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Figure  26.  Instantaneous  velocity  fields  from  BHA1R.  (a)  ^-component  and  (b)  v-component  at 
77//!  =  0.025.  Contour  interval  is  0.1;  dashed  lines  indicate  negative  values,  (c) 
Secondary  flow  at  x/X  =  0.75.  The  vectors  are  interpolated  onto  a  48x36  uniform  grid. 
The  maximum  velocity  is  0.7. 
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Figure  27.  Transverse  velocity  at  fixed  streamwise  and  vertical  position  as  a  function  of  lateral 
position  and  time:  (a)  BHA1  at  x/A  ~  0.75  and  77/A  =  0.0125 ,  contour  intervals  are  0.1. 
(b)  Case  F  at  z/H  =  0.01,  contour  intervals  are  0.05. 


Figure  29.  Transverse  velocity  variance  perturbation  from  the  wavelength  average, 
normalized  by  m*  ,  for  case  SIR.  Contour  interval  is  0.1.  Dashed  lines  indicate  negative 
values. 
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1.  INTRODUCTION 

Turbulent  dispersion  in  the  atmospheric  boundary 
layer  over  a  homogeneous  surface  is  relatively  well 
understood,  especially  under  neutral  and  convectively 
unstable  conditions.  However,  dispersion  over  terrain 
has  not  yet  been  studied  extensively.  In  this  paper, 
large-eddy  simulation  (LES),  in  conjunction  with  a 
generalized  Gaussian  puff  method,  is  used  to 
investigate  the  dispersion  of  a  passive  cloud  in  a 
convectively-driven  boundary  layer  with  moderate 
wind  speed  over  idealized  terrain.  Variations  in 
terrain  slope,  source  location  and  wind  direction  are 
considered. 

LES  is  a  modelling  technique  wherein  the  large- 
scale  eddies  are  represented  explicitly  while  small 
scale  turbulence  is  parameterized  by  a  subgrid  model. 
The  convective  boundary  layer  (CBL)  is  particularly 
amenable  to  study  by  LES,  e.g.  Mason  (1989).  LES 
has  been  used  by  Krettenaur  and  Schumann  (1992) 
and  Walko,  et  al.  (1992)  to  study  terrain-induced 
modifications  to  the  CBL,  but  with  no  mean  horizontal 
motion.  Dispersion  in  the  CBL  has  been  investigated 
using  LES,  e.g.,  Lamb  (1978)  and  Henn  and  Sykes 
(1992),  but  only  for  flat  terrain  . 

To  clearly  illustrate  the  influence  of  terrain  and 
isolate  the  effects  of  slope  and  release  location  an 
dispersion,  height  variations  are  sinusoidal  in  one 
direction  only,  obviously  an  idealized  representation 
of  ridge/valley  topography. 

2.  MODEL  DESCRIPTION 

The  LES  model  used  to  generate  the  velocity 
fields  is  an  extension  of  the  Cartesian  model  of  Sykes 
and  Henn  (1989)  and  is  based  on  the  spatially-filtered 
equations  of  motion  for  an  incompressible  Boussinesq 
fluid.  Second-order  accurate,  finite  difference  equa¬ 
tions  with  leap-frog  time  differencing  are  used  to  solve 
for  the  velocity  components  (w,v,w),  corresponding  to 
the  coordinate  axes  (x ,y,z),  and  the  perturbation 
potential  temperature.  The  subgrid  model  uses  a  tur¬ 
bulent  kinetic  energy  transport  equation  and  an  alge¬ 
braic  filter  scale.  The  Cartesian  model  equations  are 
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recast  using  the  terrain-following  coordinate  transfor¬ 
mation  of  Clark  (1977).  Thus,  the  model  vertical 
coordinate  is  defined  by 

g  =  (z-h)/J,  (1) 

where 

J  =  l-h/D.  (2) 

Here,  h(x,y)  is  the  local  terrain  elevation  and  D  is  the 
depth  of  the  computational  domain.  Details  on  the 
transformed  equations  can  be  found  in  Clark  (1977) 
and  Krettenauer  and  Schumann  (1992). 

Periodic  boundary  conditions  are  imposed  hori¬ 
zontally.  The  upper  boundary  at  z-D  is  a  rigid, 
stress-free  lid.  Rayleigh  damping  is  applied  to  the 
vertical  velocity  near  the  upper  lid  to  provide  an 
effective  non-reflecting  boundary  for  propagating 
gravity  waves.  A  uniform  heat  flux  is  specified  on  the 
bottom  surface,  z  =  h(x,y).  Monin-Obukov  similarity 
relations  are  employed  (with  suitable  rotations 
accounting  for  the  local  terrain  slope)  to  set  tempera¬ 
ture  and  velocity  on  the  first  grid  point  above  the 
surface. 

The  passive  scalar  conservation  equation  is 
solved  by  means  of  a  generalized  Gaussian  puff 
model  as  described  in  Sykes  and  Henn  (1995).  Each 
puff  has  an  associated  mass,  centroid  and  complete 
second-moment  tensor  which  gives  the  Gaussian 
spread  about  its  centroid.  The  method  accounts 
completely  for  shear  distortions  of  the  puffs  and 
employs  an  efficient  splitting/merging  scheme  to 
maintain  puff  size  consistent  with  the  LES  resolution. 
The  LES  subgrid  turbulent  kinetic  energy  is  used  to 
estimate  puff  diffusion.  A  puff  is  advected  by  the  LES 
velocity  field  interpolated  at  the  puff  centroid;  the  peri¬ 
odicity  of  the  velocity  fields  is  utilized  when  a  puff  is 
transported  beyond  the  nominal  LES  domain. 

3.  SIMULATION  CONDITIONS 

The  terrain  used  in  this  study  is  a  sinusoidal  ridge 
defined  as 

h(x,y)  =  h^(  1  -  cos2;ev//1)/2  (3) 

where  A  is  the  wavelength  and  h0  is  the  maximum 
height.  The  maximum  slope  is  given  by  S  =  7thQ/A. 
Results  are  presented  for  heights  of  0,  159m  and 
318m,  which,  with  a  wavelength  of  2km,  correspond 
to  S  =  0,  0.25  and  0.5,  respectively.  A  surface  rough 
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ness  height  of  1m  and  surface  temperature  flux  of 
O.OSKms"1  are  specified  for  all  simulations.  The 
mean  boundary  layer  depth  is  about  1000m,  giving  a 
convective  velocity  scale,  w.,  of  1ms"1. 

The  flow  is  initialized  with  a  geostrophic  wind 
(UgVg)  and  a  uniform  vertical  temperature  gradient  of 
O.OOSKm"1  above  f=800m;  a  constant  temperature 
and  horizontal  wind  are  prescribed  below. 
Turbulence  is  initiated  by  a  small  random  perturbation 
in  the  temperature  field  below  800m.  The  disper¬ 
sion  calculations  commence  only  after  an  approxi¬ 
mate  balance  is  achieved  between  the  Coriolis  force 
and  the  surface  momentum  flux.  This  balance 
typically  requires  an  integration  time  of  about  40000s. 

The  LES  calculations  use  a  48x48x65  grid  in  a 
4km  cube  domain.  The  horizontal  grid  spacing  is 
uniform.  The  vertical  grid  spacing  is  non-uniform  with 
10m  resolution  at  the  surface,  expanding  to  40m  in 
the  middle  of  the  boundary  layer. 

The  dispersion  calculations  are  initialized  from  a 
spherical  puff  10m  above  the  surface  with  a  Gaussian 
spread  of  10m.  Independent  dispersion  realizations 
are  obtained  by  varying  the  horizontal  location  of  the 
source  within  the  same  LES  fields  as  well  as  by  re¬ 
leasing  puffs  1800s  apart.  To  examine  the  effects  of 
terrain  slope,  source  location  and  wind  direction, 
seven  cases  are  considered  in  this  study  and  are 
summarized  in  Table  1.  The  dispersion  statistics  pre¬ 
sented  in  Section  5  result  from  24  realizations  of 
1800s  duration  for  each  case. 


Table  1.  Case  Summary 


Case 

5 

Source 

Location 

u> 

(ms-1) 

(ms-1) 

F 

0 

arbitrary 

5 

0 

HI 

0.25 

hill 

5 

0 

VI 

0.25 

valley 

5 

0 

H2 

0.5 

hill 

5 

0 

V2 

0.5 

valley 

5 

0 

H2D 

0.5 

hill 

0 

5 

V2D 

0.5 

valley 

0 

5 

4.  VELOCITY  FIELDS 

Since  the  primary  focus  in  this  paper  is  the 
dispersion  calculations,  only  a  brief  review  of  the 
velocity  field  results  is  presented. 

Figure  1  shows  the  mean  flow  over  ridges  with 
5  =  0.5  and  (UgVg)  =  (5,0).  The  streamlines  show  a 
clear  reversed  flow  in  the  valleys  and  acceleration 
over  the  ridge  crests,  as  expected.  The  transverse 
velocity  component,  v,  (not  shown)  is  relatively  well- 
mixed  with  little  response  to  the  local  terrain. 

An  example  of  the  instantaneous  flow  structure  is 
given  in  Figure  2.  The  near-surface  streaklines,  ob¬ 
tained  by  calculating  trajectories  in  the  instantaneous 
(«,v)-field  10m  above  the  surface,  show  the  distinct 


separation  region  between  the  two  ridges.  The  pre¬ 
dominant  flow  near  the  surface  is  in  the  positive  y- 
direction,  except  over  the  hill  crests  where  the  w- 
component  is  accelerated.  The  separation  and  reat¬ 
tachment  lines,  although  meandered  by  the  large 
eddies,  are  readily  visible  from  the  streaklines.  The 
flow  is  clearly  oriented  along  the  valley  axis  in  the 
separation  and  reattachment  regions  since  the  v- 
component  is  relatively  constant. 


Figure  1.  Mean  streamlines,  over  2km  ridges  with 
5  =  0.5.  Contour  interval  is  200m2s_1  for  positive 
values,  longs'"1  for  negative  values  (shown  dashed). 


0  12  3  4 

x  (km) 


Figure  2.  Instantaneous  streaklines  at  £=  10m  over 
2km  ridges  with  5  =  0.5. 

Vertical  profiles  of  the  three  velocity  variance 
components  are  shown  in  Figure  3  along  with  the  flat 
surface  results.  The  general  increase  with  terrain  in 
all  three  components  is  immediately  obvious.  The 
horizontal  variances  are  roughly  doubled  in  the  main 
boundary  layer,  with  significant  perturbations  near  the 
surface. 

Although  only  results  for  5  =  0.5  are  shown  here, 
those  for  the  5  =  0.25  case  are  qualitatively  similar  in 
the  sense  of  an  overall  (but  smaller)  increase  in  hori¬ 
zontal  variance.  A  notable  difference  is  the  absence 
of  mean  flow  separation  with  5  =  0.25.  (But  at  any 
given  time  there  generally  are  isolated  occurrences  of 
reverse  flow.) 


Figure  3.  Vertical  profiles  of  velocity  variances  for 
S  =  0.5  at  the  hill  crest.  The  dashed  lines  are  for  the 
corresponding  flat  case.  Both  resolved  and  total 
variances  are  shown. 


Figure  4.  Instantaneous  cloud  from  a  realization  of 
case  V2  at  840s.  (a)  Vertically-integrated  concentra¬ 
tion.  (b)  Horizontally-integrated  concentration. 
Normalized  contours  of  .01 ,  .02,  0.05,  0.1, 0.2,  0.5. 

5.  DISPERSION  STATISTICS 

A  number  of  independent  realizations  of  the 
contaminant  field  have  been  generated  with  the  goal 
of  obtaining  reliable  ensemble  dispersion  statistics. 
An  example  of  the  instantaneous  field  from  a  realiza¬ 
tion  of  case  V2  is  shown  in  Figure  4  and  clearly  illus¬ 
trates  how  terrain  can  enhance  dispersion  by  the 
separation  of  material  between  adjacent  valleys. 

To  quantify  the  effects  of  terrain  on  dispersion 
estimates  of  the  ensemble  cloud  centroid  and  spread 


for  the  cases  given  in  Table  1  will  now  be  presented. 
The  cloud  centroid  for  each  realization  is  defined  as 

where  the  summation  is  over  all  puffs,  Qn  is  the  mass 
of  puff  n  and  xln  is  the  puff  centroid  relative  to  the  re¬ 
lease  location.  The  ensemble  centroid  is  then  {x^ 
where  the  brackets  indicate  averaging  over  all  real¬ 
izations.  The  total  ensemble  spread,  Taa ,  is  the  com¬ 
bination  of  a  meander  component,  Maa  and  a  relative 
dispersion  component,  Saa : 

Taa  =  Maa+Saa  (5) 

(a=  1,  2  or  3  with  no  summation,  corresponding  to 
coordinates  x,  y,  or  z)  where 


and 


(7) 


Here,craan  are  the  diagonal  components  of  the  sec¬ 
ond-moment  tensor  for  each  puff. 

The  mean  horizontal  trajectories  of  the  ensemble 
clouds,  relative  to  the  release  locations,  are  shown  in 
Figure  5.  For  the  Ug=  5ms“1  cases,  the  paths 
become  increasingly  oriented  parallel  to  the  ridges 
with  increasing  slope,  a  result  of  the  Coriolis  force 
balancing  the  increased  surface  stress  normal  to  the 
ridges.  The  valley  release  with  S  =  0.5,  case  V2,  is 
particularly  interesting  since  its  path  initially  shows  a 
slight  reverse  flow  in  the  x-direction  and,  as  a  result,  is 
confined  to  the  valley  and  travels  along  its  axis  for 
roughly  the  first  500s,  after  which  it  parallels  the  hill 
release.  On  the  other  hand,  for  S  =  0.25,  the  paths  for 
the  hill  and  valley  releases  diverge  only  slightly  at  the 
outset  and  then  quickly  converge. 


Figure  5.  Ensemble  cloud  centroid  trajectories  for  the 
cases  of  Table  1.  Symbols  represent  120s  time  in¬ 
crements.  The  flat  case  F*  is  shown  for  comparison 
with  the  Vg  =  Sms"1  cases  and  is  merely  case  F 
rotated  90°. 

The  ensemble  cloud  spread  in  the  direction  nor¬ 
mal  to  the  ridges  as  a  function  of  time  is  shown  in 


(10) 


Figure  6  for  the  first  5  cases  from  Table  1,  i.e.,  Ug  = 
5ms'1.  It  can  be  seen  that  Txx  is  considerably  en¬ 
hanced  for  S  =  0.5  over  the  flat  terrain  case  for  both 
hill  and  valley  releases.  And  while  the  hill  release  with 
5  =  0.25  also  shows  enhanced  dispersion,  the  valley 
release  surprisingly  results  in  reduced  dispersion. 


time  (s) 

Figure  6.  Ensemble  dispersion  normal  to  the  ridges 
for  the  t/^  =  5ms'1  cases.  Solid  line:  case  F;  short 
dashes:  HI;  long  dash,  short  dash:  VI;  long  dashes: 
H2;  long  dash,  two  short  dashes:  V2.  (a)  Total,  (b) 
relative  dispersion,  (c)  meander  component. 

Sensitivity  to  release  location  for  the  terrain 
cases  is  not  surprising.  V2  shows  greater  dispersion 
normal  to  the  ridges  than  H2  and  this  is  most  likely 
related  to  the  flow  separation  in  the  valleys  generating 
a  significantly  larger  meander  component  for  V2  than 
for  H2.  Figure  5b  clearly  shows  for  V2  and  H2 
diverging  over  the  first  600s,  after  which  they  slowly 
converge.  Thus,  sensitivity  to  release  location  is 
mainly  due  to  the  different  time  and  spatial  scales  of 
the  large  eddies  at  the  hill  crest  and  in  the  valley. 

The  dispersion  for  case  VI  appears  to  be 
anomalous  and  is  clearly  due  to  the  behavior  of 
since  5^  is  quite  close  to  cases  F  and  HI.  After 
closely  matching  the  V2  results  for  the  first  400s,  the 
VI  meander  component  decreases  slightly  and  then 
remains  nearly  constant  for  the  remainder  of  the  inte¬ 
gration,  showing  no  apparent  convergence  with  the 
corresponding  hill  release.  VI  is  the  only  case  which 
shows  any  decrease  or  constancy  of  M^;  all  others 
show  monotonic  increases.  The  lack  of  convergence 
with  HI  is  especially  surprising  given  that  the  mass- 
weighted  (resolved)  velocity  variance,  Tuu ,  does  con¬ 
verge  for  the  two  cases,  as  shown  in  Figure  7.  Here 
Tuu  =  Muu  +  Suu  (8) 

where 

mu„=((«-{“»2)  (9) 


un  is  the  w-component  at  the  centroid  of  puff  n  and  u 
is  the  mass-weighted  mean  for  a  realization.  It  can  be 
seen  that  Tuu  for  VI  is  initially  similar  to  V2,  then 
decreases  at  a  time  corresponding  to  the  decrease  in 
Tja  and  finally  converges  with  case  HI  at  around 
800s.  The  fact  that  does  not  converge  for  VI  and 
HI  despite  this  implies  some  kind  of  negative  velocity 
correlation  for  the  valley  release.  Since 


dM„ 

dt 


=  2  (nr) 


(ii) 


where  the  prime  denotes  fluctuation  from  the  ensem¬ 
ble  mean,  Mxv  can  decrease  or  remain  constant  only  if 
(FF)  <  0  and  indeed  this  condition  occurs  for  VI  but 
not  the  other  cases.  It  is  unknown  if  this  is  a  general 
result  or  an  artifact  of  the  sinusoidal  character  of  the 
terrain  and  clearly  needs  further  investigation. 


fluctuation  intensities  for  the  Ug  =  5ms  1  cases.  Line 
patterns  as  in  Figure  6. 


time  (s) 

Figure  8.  Ensemble  dispersion  parallel  to  the  ridges 
for  the  Ug  =  5ms'1  cases.  Line  patterns  as  in  Figure  6. 


Dispersion  parallel  to  the  ridges  is  shown  in 
Figure  8.  All  terrain  cases  show  enhanced  dispersion 
over  the  flat  case  but,  with  the  exception  of  the 
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anomalous  case  VI ,  the  absolute  values  are  less  than 
in  the  normal  direction.  The  enhanced  dispersion  is 
mostly  due  to  larger  meander  components;  relative 
dispersion  in  the  presence  of  terrain  is  not 
significantly  greater  than  the  flat  case.  Also,  as 
expected,  sensitivity  to  release  location  is  not  as  pro¬ 
nounced  as  in  the  ^-direction  since  v  does  not  vary 
significantly  with  the  local  terrain. 


0  400  800  1200  1600 

Figure  9.  Ensemble  dispersion  normal  to  the  ridges 
for  the  Vg  =  Sms-1  cases.  Solid  line:  case  F  (rotated); 
long  dashes:  H2D;  long  dash,  two  short  dashes:  V2D. 
(a)  Total,  (b)  relative  dispersion,  (c)  meander 
component. 


Figure  10.  Trajectories  of  cloud  centroids  from  indi¬ 
vidual  realizations  of  cases  V2D  and  H2D. 

Finally,  to  show  that  the  orientation  between  the 
wind  and  predominant  terrain  direction  can  be  impor¬ 
tant,  statistics  for  cases  V2D  and  H2D  are  now  pre¬ 
sented.  Figure  5  shows  the  horizontal  ensemble  tra¬ 
jectories  of  these  cases  in  addition  to  the  flat  case  F, 
rotated  90  for  comparison.  Obviously,  the  valley 
release  results  in  slower  transport  across  the  ridges, 
implying  that  material  may  initially  be  trapped  within  a 


valley.  However,  as  shown  in  Figure  9,  the  dispersion 
across  the  ridges  is  greater  for  V2D  than  H2D  due  to 
an  increased  meander  component  for  the  valley 
release.  Apparently,  the  cloud  in  some  realizations  of 
V2D  remains  in  the  valley  in  which  it  was  released  for 
a  relatively  long  time,  while  in  others,  the  cloud 
quickly  moves  over  the  ridge.  Examples  of  this 
behavior  are  given  in  Figure  10  which  shows  the  cen¬ 
troid  trajectories  of  all  realizations  of  cases  V2D  and 
H2D. 

6.  CONCLUSION 

LES  has  been  used  to  investigate  the  flow  and 
dispersion  over  idealized  ridge/valley  terrain,  it  has 
been  found  that  dispersion  is  generally  enhanced  by 
terrain,  especially  for  large  slopes,  and  that  it  can  be 
sensitive  to  release  location.  It  is  also  found  that  dis¬ 
persion  is  a  function  of  wind  direction  with  respect  to 
the  ridge/valley  orientation,  greater  dispersion  result¬ 
ing  from  flow  predominantly  normal  to  the  ridge.  An 
anomalous  decrease  in  dispersion  (relative  to  flat  ter¬ 
rain)  for  a  valley  release  and  moderate  slope  certainly 
requires  further  investigation. 

Future  work  will  be  aimed  at  quantifying  the  dis¬ 
persion  results  in  terms  of  terrain  slope,  wind  speed 
and  orientation  and  extending  this  to  more  complex 
terrain. 
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1.  INTRODUCTION 

Most  of  our  current  understanding  of  dispersion 
in  the  atmosphere  is  based  on  an  assumption  of  hori¬ 
zontal  surface  homogeneity,  which  is  valid  in  only 
limited  situations.  Since  significant  terrain  features 
cover  much  of  the  earth’s  surface,  dispersion  in  the 
presence  of  terrain  is  of  considerable  interest.  Here, 
we  use  large-eddy  simulation  (LES)  to  investigate  the 
dispersion  of  a  passive  cloud  over  idealized  terrain. 
Variations  in  terrain  slope  and  source  location  under 
strongly  convective  and  near-neutral  conditions  are 
considered. 

LES  is  a  modeling  technique  wherein  the  large- 
scale  eddies  are  represented  explicitly  while  a  subgrid 
model  parameterizes  small-scale  turbulence.  The 
convective  boundary  layer,  characterized  by  eddies 
on  the  scale  of  the  boundary  layer,  is  particularly 
amenable  to  study  by  LES,  e.g.  Mason  (1989).  Like¬ 
wise,  the  presence  of  large-scale  terrain  features, 
while  certainly  presenting  modeling  difficulties,  is  ex¬ 
pected  to  generate  flow  features  well-represented  by 
LES  (Gong,  et  al.,  1996).  Terrain-induced  modifica¬ 
tions  to  the  convective  boundary  layer  have  been 
studied  using  LES  by  Krettenaur  and  Schumann 
(1992),  and  Dornbrack  and  Schumann  (1993). 

To  clearly  illustrate  the  influence  of  terrain  and 
isolate  the  effects  of  slope  and  release  location  on 
dispersion,  height  variations  are  sinusoidal  in  one 
direction  only.  Although  obviously  an  idealization, 
this  is  a  reasonable  prototype  of  ridge-valley  terrain. 

2.  MODEL  DESCRIPTION 

The  LES  model  used  to  generate  the  velocity  fields  is 
an  extension  of  the  Cartesian  model  of  Sykes  and 
Henn  (1989)  and  is  based  on  the  spatially  filtered 
equations  of  motion  for  an  incompressible  Boussi- 
nesq  fluid.  Second-order  accurate,  finite  difference 
equations  with  leap-frog  time  differencing  are  used  to 
solve  for  the  resolved-scale  velocity  components 
(u,v,w ),  corresponding  to  the  coordinate  axes  (x,y,z), 
and  the  perturbation  potential  temperature.  The  sub- 

*  Corresponding  author  address :  D.  S.  Henn,  Titan 
Corp.,  ARAP  Group,  50  Washington  Road,  Princeton, 
NJ  08543-2229;  e-mail:  dhenn@titan.com 


grid  model  uses  a  turbulent  kinetic  energy  transport 
equation  and  an  algebraic  filter  scale.  The  Cartesian 
model  equations  are  recast  using  a  terrain-following 
coordinate  transformation;  details  can  be  found  in 
Clark  (1977)  and  Krettenauer  and  Schumann  (1992). 
Preliminary  results  from  this  model  have  been  re¬ 
ported  in  Henn  and  Sykes  (1996,  1997). 

Periodic  boundary  conditions  are  imposed  hori¬ 
zontally.  The  upper  boundary  is  a  rigid,  stress-free 
lid.  Rayleigh  damping  is  applied  to  the  vertical  veloc¬ 
ity  near  the  upper  lid  to  provide  an  effectively  non¬ 
reflecting  boundary  for  propagating  gravity  waves.  A 
uniform  heat  flux  is  specified  on  the  bottom  surface. 
Monin-Obukov  similarity  relations  are  employed  (with 
suitable  rotations  accounting  for  the  local  terrain 
slope)  to  set  temperature  and  velocity  on  the  first  grid 
point  above  the  surface. 

The  scalar  dispersion  calculations  employ  the 
generalized  Gaussian  puff  model  of  Sykes  and  Henn 
(1995)  when  early-time,  near-source  behavior  is  of 
primary  interest  or  a  particle  method,  which  accounts 
for  subgrid  diffusivity  with  a  random-walk  model, 
when  late-time  statistics  are  desired.  The  second- 
order  accurate  Adams-Bashforth  scheme  is  used  in 
both  methods  to  advance  the  puff/particle  location, 
with  the  advection  velocity  determined  by  interpola¬ 
tion  from  the  LES  velocity  field.  The  periodicity  of  the 
velocity  fields  is  utilized  when  a  puff  or  particle  is 
transported  beyond  the  nominal  LES  domain. 

3.  SIMULATION  CONDITIONS 

The  terrain  used  in  this  study  is  a  sinusoidal  ridge 
defined  as 

h(xyy)-hQ(l-cos27Dc/X)/2  (3) 

where  A  is  the  wavelength  and  h0  is  the  maximum 
height.  The  maximum  slope  is  given  by  S  =  7di0/A. 
Results  are  presented  for  heights  of  0,  159  m  and 
318  m,  which,  with  A  =  2  km  correspond  to  S  =  0,  0.25 
and  0.5,  respectively.  A  surface  roughness  height  of 
1  m  is  assumed  for  all  simulations. 

The  flow  is  driven  by  a  geostrophic  wind  of 
5  ms"1  perpendicular  to  the  hill  crests.  To  investigate 
the  effects  of  thermal  stability,  two  uniform  surface 
temperature  fluxes  are  considered:  0.03°Kms"1  and 
0.0045°Kms"1.  The  resulting  convective  boundary 
layer  is  capped  by  an  overlying  stable  layer  with  a 
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mean  temperature  gradient  of  0.005°Knrn1.  The 
boundary  layer  depth  is  maintained  at  approximately 
1000  m,  which  means  that  the  two  heat  fluxes  con¬ 
sidered  correspond  to  convective  velocity  scales  of 
1  ms-1  and  0.5  ms"1,  respectively.  The  higher  heat 
flux  results  in  flow  which,  at  least  over  flat  terrain,  is 
dominated  by  large-scale  convective  eddies  (the  ratio 
of  the  shear-stress  velocity  to  the  convective  velocity 
scale, um/w.  =0.4)  while  the  lower  heat  flux  case  is 
closer  to  neutral  stability  in  character  ~  0.7). 

The  LES  calculations  use  a  48x48x65  grid  in  a 
4  km  cube  domain.  The  horizontal  grid  spacing  is 
uniform.  The  vertical  grid  spacing  is  non-uniform  with 
10  m  resolution  at  the  surface,  expanding  to  40  m  in 
the  middle  of  the  boundary  layer. 

The  short-range  (puff)  dispersion  calculations  are 
initialized  from  a  spherical  puff  10  m  above  the  sur¬ 
face  with  a  Gaussian  spread  of  10  m.  The  late-time 
(particle)  calculations  are  initialized  from  a  cloud  of 
1000  randomly  distributed  particles  located  200  m 
above  the  wave  crest  with  a  spread  of  100  m.  Inde¬ 
pendent  dispersion  realizations  within  the  same  LES 
field  are  obtained  by  separating  four  releases  1000  m 
apart;  this  procedure  is  repeated  for  a  number  of  in¬ 
dependent  LES  fields.  Thus,  the  dispersion  statistics 
presented  in  Section  5  are  computed  from  ensembles 
of  12  to  24  realizations. 

To  examine  the  effects  of  terrain  slope,  source 
location  and  stability,  sixteen  cases  are  considered 
and  are  summarized  in  Table  1.  The  case  designa¬ 
tions  are  as  follows:  “A”  and  “B”  denote  low  heat  flux 
or  high  heat  flux  cases,  respectively;  0,  1  and  2  corre¬ 
spond  to  S  -  0,  0.25  and  0.5,  respectively;  “H”  and  “V” 
denote  hill  and  valley  releases,  respectively,  while  “L” 
denotes  a  late-time  calculation. 


Table  1.  Case  Summary. 
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4.  VELOCITY  FIELDS 


Since  the  primary  focus  of  this  paper  is  the  ef¬ 
fects  of  significant  terrain  on  dispersion,  only  a  brief 
review  of  the  velocity  field  statistics  will  be  presented. 

Figure  1  shows  the  mean  streamlines  with 
5  =  0.5  for  the  two  heat  fluxes  (the  A2  and  B2  cases). 
The  streamlines  show  a  clear  reversed  flow  in  the 
valleys  and  acceleration  over  the  ridge  crests. 
Clearly,  the  separation  zone  is  smaller  for  the  higher 
heat  flux  case,  most  likely  a  result  of  greater  vertical 


mixing.  The  zero  streamline  is  an  indication  of  where, 
to  some  extent,  the  flow  in  the  valleys  is  de-coupled 
from  the  flow  above.  This  naturally  has  significant 
implications  for  dispersion.  For  example,  material 
released  in  a  valley  is  more  likely  to  become  trapped 
in  a  recirculating  flow  under  lower  heat  flux  condi¬ 
tions.  Furthermore,  less  vertical  mixing  implies  that 
the  material  will  remain  in  the  recirculation  zone  for  a 
longer  time. 


Figure  1.  Mean  streamlines  over  5  =  0.5  ridges 
with  (a)  wm  =  0.5  ms'1 ,  (b)  =  1  ms"1.  Contour  inter¬ 

val  is  40  m2s-1  for  positive  values,  10  m2s_1  for 
negative  values  (shown  dashed). 


The  5  =  0.25  cases  (not  shown)  have  much  less 
mean  flow  separation,  particularly  with  the  higher 
heat  flux.  However,  reversed  flow  still  occurs  inter¬ 
mittently. 

Vertical  profiles  of  the  mean  streamwise  velocity 
component  are  shown  in  Figure  2  for  5  =  0.5  and 
5  =  0.  The  profiles  are  presented  as  perturbations 
from  the  mean  sectional  bulk  average,  defined  as 


1 


H-h(x) 


a 

judz 


h(x) 


where  H  is  the  boundary  layer  depth.  The  profiles 
over  the  valley  are  typical  of  much  of  the  flow  and 
show  significantly  enhanced  shear  compared  with  the 
flat  cases.  The  shear  is  generally  greater  in  magni¬ 
tude  as  well  as  vertical  extent.  This  is  significant 
since  Taylor’s  analysis  for  the  longitudinal  dispersion 
coefficient,  D„,  is  based  on  a  convolution  integral  of 
the  velocity  perturbations,  namely 


D„  =  ~jujK7'ju'd&Vdz 

n  n  n 


0  0  0 

where  Kz  is  the  vertical  diffusivity  (Elder  1959). 
Hence,  the  greater  the  shear,  the  greater  the  longitu¬ 
dinal  dispersion.  Furthermore,  greater  vertical  mixing 
reduces  longitudinal  dispersion.  Although  Taylor’s 
analysis  is  only  valid  at  late  times  (after  complete 
vertical  mixing)  and  over  homogenous  terrain,  we  still 
expect  the  enhanced  shear  and  reduced  vertical 
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Figure  2.  Mean  perturbation  velocity  profiles.  Solid 
lines  indicate  flat  terrain,  dashed  lines  S  =  0.5;  thin 
lines  indicate  w*  =  0.5  ms"1,  thick  lines  w*  =  1  ms'1 . 
(a)  Hill  crest  and  (b)  valley  bottom. 

mixing  of  the  lower  heat  flux  case  to  result  in  en¬ 
hanced  dispersion. 

5.  DISPERSION  STATISTICS 

Ensemble  statistics  from  the  puffs  or  particles  were 
generated  at  selected  times  (at  intervals  of  120s  for 
the  puffs  and  100s  for  the  particles).  Due  to  space 
limitations,  only  longitudinal  statistics  (in  the  x- 
direction,  normal  to  the  hill  crests)  will  be  presented. 

The  total  ensemble  longitudinal  second  moment, 
7^,  is  the  sum  of  a  meander  component,  and  a 
relative  dispersion  component,  Sxx: 

T xx  ~  Mxx  +  $ xx 

where 

and 

Sx,=(<ria+{x-xf^ 

where  an  overbar  denotes  a  mass-weighted  average 
over  all  puffs  for  a  particular  realization  and  angle 
brackets  denotes  averaging  over  all  realizations. 
Thus,  the  ensemble  longitudinal  cloud  centroid  is 
given  by  (I).  Note  that  for  particle  statistics  the  over¬ 
bar  simply  indicates  averaging  over  all  particles  of  a 
realization  since  no  mass  is  associated  with  a  parti¬ 
cle.  Here,  is  the  individual  puff  longitudinal  sec¬ 
ond-moment  about  its  centroid  (zero  for  a  particle). 

These  moments,  as  functions  of  time,  are  shown 
in  Figure  3  for  the  various  short-range  cases.  It  can 
be  seen  that  the  total  spread  7^  is  quite  similar  for  A0 
and  B0,  but  this  is  actually  due  to  offsetting  trends  in 
the  relative  and  meander  components.  Not  surpris¬ 
ingly,  is  greater  in  B0  due  to  the  large-scale  ed¬ 
dies  characteristic  of  the  convective  boundary  layer. 
The  larger  S ^  in  AO  is  probably  due  to  the  reduced 
vertical  mixing,  which  allows  puffs  to  remain  in  the 
near-surface  shear  region  for  a  relatively  long  time. 


The  effects  of  terrain,  at  least  for  the  slopes  con¬ 
sidered  in  this  study,  are  quite  substantial  and  show 
sensitivity  to  release  location  and  thermal  stability. 
The  S  =  0.25  hill  releases  (A1H  and  B1H)  probably 
show  the  smallest  differences  from  the  flat  terrain 
cases,  but  nonetheless  some  subtle  effects  are  ap¬ 
parent.  For  approximately  the  first  500  s,  the  total 
spread  is  reduced  compared  to  A0  and  B0  due  to  a 
reduction  in  the  meander  component.  (This  is  also 
seen  in  A2H  and  B2H.)  Subsequently,  the  relative 
component  in  A1H  is  enhanced  relative  to  both  A0 
and  B1 H,  most  likely  due  to  terrain-induced  shear. 

The  valley  releases  with  S  =  0.25  show  rather 
surprising  behavior  regarding  the  meander  compo¬ 
nent.  In  both  A1V  and  B1V,  increases  rather 
rapidly  for  the  first  500  s  and  then  levels  off.  The 
early  behavior  is  probably  due  to  the  intermittency  of 
the  separation  in  the  valleys:  while  some  releases 
initially  get  trapped  in  a  recirculation  zone  and  move 
slightly  upstream  for  some  time,  others  simply  move 
downstream,  thus  resulting  in  a  wide  separation  be¬ 
tween  centroids  of  different  realizations.  The  fact  that 
Mm  then  levels  off  implies  that  the  centroids  of  differ¬ 
ent  realizations  move  downstream  more  or  less  in 
parallel.  This,  in  turn,  is  probably  an  indication  that 
little  material  is  becoming  trapped  in  a  neighboring 
valley  (which  we  assume  would  occur  intermittently 
and  hence  result  in  an  increase  in  M^).  In  contrast, 
puffs  from  some  of  the  hill  releases  do  become 
trapped,  especially  since  their  small  size  at  this  stage 
means  they  can  easily  be  transported  by  downdrafts 
and  recirculation  zones.  Thus,  especially  in  B1H,  the 
meander  component  is  still  increasing  at  1800  s.  This 
may  also  partially  explain  the  larger  relative  compo¬ 
nents  for  A1 H  and  B1 H:  will  likely  increase  if  some 
material  gets  trapped  in  a  separation  bubble. 

All  S  =  0.5  cases  ultimately  result  in  enhanced 
total  dispersion  compared  with  the  corresponding 
S  =  0.25  and  flat  cases.  This  is  principally  due  to  the 
relative  components  which  are  enhanced  both  by  the 
increase  in  mean  shear  and  possibly  by  material  get¬ 
ting  trapped  in  recirculating  flow.  As  with  the  S  =  0.25 
cases,  is  apparently  converging  for  the  hill  and 
valley  releases  by  the  end  of  the  calculations.  Early 
on,  Sxx  is  somewhat  larger  for  the  hill  releases,  as  with 
A1H  and  B1H,  due  to  material  being  brought  down 
immediately  into  a  valley  and  getting  sheared  out. 

The  meander  components  of  A2V  and  B2V  are 
an  interesting  study  in  contrasts.  in  B2V  behaves 
in  a  manner  similar  to  A1V  and  B1V  in  that  it  in¬ 
creases  rapidly  for  approximately  the  first  500  s. 
Even  with  the  larger  slope,  the  separation  in  this  case 
is  rather  intermittent,  which  as  explained  above,  re¬ 
sults  in  a  relatively  large  meander  component.  How¬ 
ever,  after  the  first  200  s,  in  A2V  is  reduced  rela¬ 
tive  to  A2H,  B2H  and  B2V.  In  this  case,  the  separa¬ 
tion  zones  are  larger  and  tend  to  persist  longer  than 
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Figure  3.  Longitudinal  spreads  from  puff  calculations  for  various  hill  slopes  and  release  locations.  Solid  lines 
indicate  hill  releases,  dashed  lines  valley  releases;  thin  lines  correspond  to  "A"  cases,  thick  lines  to  "B"  cases. 


for  the  lower  heat  flux  case.  Thus  virtually  all  realiza¬ 
tions  of  A2V  are  initially  released  into  a  recirculation 
zone  and  remain  in  it  for  some  time,  with  the  result 
that  the  separation  of  the  individual  realizations  is 
reduced.  This  is  illustrated  in  Figure  4,  which  shows 
the  longitudinal  centroid  locations  for  individual  reali¬ 
zations  of  A2V  and  B2V.  Note  that  the  centroids  of  all 
A2V  realizations  remain  in  the  original  valley 
(I<1  km)  for  the  first  1000  s;  about  half  are  still  in  at 
1800  s,  the  end  of  the  calculation.  In  contrast,  while 
most  realizations  of  B2V  do  evidently  get  released 
into  recirculating  flow,  many  nonetheless  leave  the 
valley  quite  quickly  and  all  have  left  by  approximately 
1400  s.  Figure  4  also  illustrates  why  in  A2V  is 
still  increasing  at  the  end  of  the  calculation  while  it  is 
leveling  off  in  B2V.  The  A2V  realizations  show  a 
large  variation  in  valley  residence  times,  resulting  in  a 
rapid  increase  of  the  meander  component.  This  will 
continue  until  all  realizations  have  left  the  valley,  as 
has  already  occurred  with  B2V. 

We  now  examine  longer  simulations  where  the 
clouds  become  quite  large,  typically  spanning  a  num¬ 
ber  of  hill  wavelengths,  and  presumably  the  initial 
source  conditions  are  unimportant.  Since  a  puff  cal¬ 
culation  would  require  a  prohibitive  number  of  puffs, 
we  use  particles  to  generate  these  dispersion  statis¬ 
tics.  While  the  particle  simulations  cannot  provide 
detailed  information  about  the  concentration  fields, 
they  are  adequate  for  computing  centroids  and  sec¬ 
ond-moments.  We  are  mainly  interested  in  deter¬ 
mining  “late-time”  behavior,  which  generally  implies 


steady  growth  of  the  second-moment  and  can  only 
exist  when  material  is  vertically  well-mixed. 

The  longitudinal  second-moments  from  these 
cases  are  shown  Figure  5.  Generally,  the  spreads 
increase  with  slope.  This  is  due  entirely  to  the  rela¬ 
tive  component  since,  as  expected,  the  meander 
components  all  level  off  after  a  while.  The  S  =  0.5 
cases  A2L  and  B2L  take  the  longest  to  obtain  a  con¬ 
stant  Mja,  most  likely  due  to  the  larger-scaie  motions 
induced  by  the  larger  terrain  features. 

The  increase  in  S ^  with  slope  is  most  evident  in 
the  lower  heat  flux  cases.  In  contrast,  the  higher  heat 
flux,  moderate  slope  case,  B1L,  shows  virtually  no 
enhanced  dispersion  over  B0L,  although  B2L  does 
show  considerably  larger  relative  dispersion.  For 
both  non-zero  slopes  considered,  the  lower  heat  flux 


time  (s) 

Figure  4.  Longitudinal  centroid  position  of  individual 
realizations  from  (a)  A2V  and  (b)  B2V. 
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Figure  5.  Longitudinal  spread  from  late-time  particle 
simulations.  Solid  lines:  AOL  (thin),  BOL  (thick);  long 
dashes:  AIL  (thin),  B1L  (thick);.  Short  dashes:  A2L 
(thin),  B2L  (thick). 


Figure  6.  Longitudinal  diffusion  coefficients  from  late¬ 
time  calculations.  Line  patterns  as  in  Figure  5. 


cases  show  greater  dispersion  than  the  correspond¬ 
ing  high  heat  flux  cases. 

It  is  convenient  to  define  a  time-dependent  lon¬ 
gitudinal  dispersion  coefficient  based  on  the  rate  of 
increase  in  the  longitudinal  second  moment: 


This  becomes  equivalent  to  Taylor’s  late-time  disper¬ 
sion  coefficient  when  dS^/dt  becomes  constant. 
(We  use  the  relative  component  since  Figure  5  shows 
that  the  meander  component  is  negligible  at  late¬ 
time.)  Figure  6  shows  that  D  is  greatly  enhanced  by 
the  terrain  in  this  study,  except  for  B1L,  which  is  not 
very  different  from  the  flat  terrain  case.  Most  cases 
have  achieved  more  or  less  steady  values  by  ap¬ 
proximated  5000  s,  the  exception  being  B2L.  A2L 


shows  quite  rapid  early  growth  before  leveling  off  to 
D„  =  4000  mV1 ,  roughly  eight  times  the  flat  terrain 
value.  The  continued  growth  in  B2L  is  rather  sur¬ 
prising  given  the  length  of  this  integration.  It  may 
reflect  some  unsteadiness  in  the  wind  fields  on  the 
Coriolis  timescale  and  needs  further  investigation. 

6.  CONCLUSION 

LES  has  been  used  to  investigate  dispersion 
over  idealized  ridge/valley  terrain.  It  has  been  found 
that  dispersion  is  generally  enhanced  by  terrain,  but 
that  it  can  be  sensitive  to  release  location  with  respect 
to  flow  features  such  re-circulation  zones,  intermittent 
separation,  vertical  mixing  and  velocity  shears.  Late¬ 
time  dispersion  is  shown  to  be  significantly  enhanced 
for  large  slopes,  the  dispersion  coefficient  being  eight 
times  the  flat  terrain  value  for  the  largest  slope  and 
lowest  heat  flux  considered  in  this  study. 
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